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Abstract 


Unsteady  flow  separation  during  dynamic  stall  often  leads  to  unacceptably  large  vibra¬ 
tory  loads  and  acoustic  noise,  and  limit  forward  flight  speeds  and  maneuverability.  To 
gain  a  quantitative  understanding  of  the  unsteady  separation  process,  large-eddy  simula¬ 
tion  (LES)  of  turbulent  flow  over  a  pitching  airfoil  at  realistic  Reynolds  and  Mach  numbers 
is  performed.  Numerical  stability  at  high  Reynolds  number  simulation  is  maintained  us¬ 
ing  an  unstructured-grid  LES  technology,  which  obeys  higher-order  conservation  principles 
and  employs  a  global-coefficient  subgrid-scale  turbulence  model.  A  hybrid  implicit-explicit 
time-integration  scheme  is  employed  to  provide  a  highly  efficient  way  to  treat  time-step  size 
restriction  in  the  separated  flow  region  locally  refined  with  dense  mesh.  The  present  simula¬ 
tions  confirm  the  stability  and  effectiveness  of  the  presented  numerical  schemes  for  dynamic 
stall  simulations  at  realistic  operating  Reynolds  and  Mach  numbers  and  show  the  charac¬ 
teristics  of  flow  separation  and  reattachment  processes  which  are  qualitatively  congruent 
with  experimental  observation. 

To  improve  quantitative  understanding  of  unsteady  separation  processes  of  turbulent 
boundary  layers,  direct  numerical  simulations  are  performed.  The  distinct  characteristics 
of  unsteady  separating  turbulent  boundary  layers  are  revealed  by  a  systematic  comparison 
with  steady  attached/separated  turbulent  boundary  layers.  For  this  purpose,  four  different 
flow  configurations  are  simulated:  (i)  turbulent  boundary  layer  flow  with  a  zero-pressure 
gradient;  (ii)  turbulent  boundary  layer  flow  with  an  adverse-pressure  gradient;  (iii)  steady 
separated  turbulent  boundary  layer  flow;  and  (iv)  unsteady  separating  turbulent  boundary 
layer  flow.  The  present  comparative  study  suggests  physical  phenomena  during  the  unsteady 
separation  process  including  unsteady  boundary-layer  detachment  and  reattachment,  and 
production  and  dissipation  of  turbulent  kinetic  energy  and  vorticity. 
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Chapter  1 


Introduction 

Dynamic  stall  is  a  nonlinear  and  unsteady  aerodynamic  phenomenon  resulting  in  stall  delay 
during  a  time-dependent  motion  of  an  airfoil  at  angles  of  attack  higher  than  its  static  stall 
angle.  Dynamic  stall  of  high  Army  relevance  occurs  on  the  retreating  blade  of  a  helicopter 
rotor  experiencing  a  pitching  motion  which  leads  to  unsteady  flow  separation  followed  by 
load  and  pitching-moment  overshoots.  The  unsteady  flow  separation  can  in  turn  lead  to 
unacceptably  large  vibratory  loads  and  acoustic  noise,  and  limit  forward  flight  speeds,  load, 
and  maneuverability  (Lorber  &  Carta  1994).  The  unsteady  separation  is  reported  to  be 
influenced  by  the  Reynolds  and  Mach  numbers,  blade  geometry,  pitch  rate,  and  freestream 
turbulence  level  (McCroskey  1982). 

Numerous  investigations  of  unsteady  separation  associated  with  dynamic  stall  have  been 
conducted  at  chord-based  Reynolds  numbers  ( Rec )  in  the  range  of  103  —  10',  at  Mach 
numbers  for  incompressible  to  transonic  flow,  and  for  a  wide  variety  of  blade  geometries. 
Most  experimental  studies  have  concentrated  on  measurements  of  aerodynamic  forces  such 
as  the  surface  pressure  and  overall  loads  ( e.g .,  Lorber  &  Carta  (1987);  Jumper  et  al.  (1987)), 
or  on  the  flow  field  visualization  (e.g.,  McCroskey  et  al.  (1976)).  Quantitative  measurements 
of  the  separated  flow  field  and  wake  around  a  pitching  airfoil  have  been  difficult  using 
experimental  techniques,  and  therefore,  have  rarely  been  reported  in  the  literature. 

Computational  fluid  dynamics  has  become  increasingly  useful  in  studying  dynamic  stall 
(see  Ekaterinaris  &;  Platzer  (1997)  for  a  review).  Computational  works  have  often  been 
performed,  especially  at  practical  Reynolds  numbers,  using  the  Reynolds- averaged  Navier- 
Stokes  (RANS)  equations  or  its  unsteady  counterpart  (URANS)  (e.g.,  Visbal  (1988);  Spent- 
zos  et  al.  (2005)).  However,  it  is  known  to  be  very  challenging  for  (U)RANS  to  accurately 
predict  highly  unsteady  flow  involving  incipient  flow  separation,  formation  and  evolution 
of  stall  vortices,  and  reattachment.  So  far,  none  of  the  available  turbulence  models  has 
shown  a  satisfactory  predictive  performance  for  all  of  these  flow  phenomena.  In  a  review 
article,  Carr  &  McCroskey  (1992)  concluded  that  “Turbulence  modeling  becomes  of  crucial 
importance  when  dynamic  stall  is  considered.  This  is  particularly  true  when  the  question  of 
incipient  separation  and  dynamic-stall-vortex  development  is  to  be  represented  by  a  single 
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turbulence  model;  under  this  condition,  the  use  of  a  turbulence  model  based  on  equilibrium 
attached  boundary  layers  in  steady  flow  (e.g.,  eddy  viscosity,  Baldwin-Lomax)  is  open  to 
question.  The  task  of  predicting  separation  by  definition  deals  with  boundary  layers  that 
have  experienced  very  strong  pressure  gradients,  often  both  positive  and  negative;  the  flow 
approaching  unsteady  separation  contains  high  levels  of  vorticity  induced  by  these  pres¬ 
sure  gradients,  and  is  strongly  unsteady.  Recent  study  has  shown  that  modification  of  the 
turbulence  model  can  completely  change  the  resulting  flow  results;  at  the  same  time,  very 
little  has  been  experimentally  documented  about  the  character  of  turbulence  under  these 
conditions.” 

The  detached-eddy  simulation  (DES)  and  kinetic-energy  simulation  (KES,  Fang  Sz  Menon 
2006)  have  also  been  explored  in  dynamic  stall  applications.  Although  encouraging  results 
were  obtained  using  DES  over  a  wide  range  of  flow  configurations,  flow  separation  is  often 
difficult  to  predict  due  to  the  existence  of  the  gray  zone  in  which  the  model  is  improperly 
converted  from  attached  boundary  layer  to  massive  separation  (Spalart  et  al.  1997).  KES, 
which  was  designed  to  predict  flow  at  length  scales  between  the  computational  grid  scale 
and  the  integral  scale,  is  in  the  evaluation  stage  as  discussed  by  its  developers  (Fang  Sz 
Menon  2006). 

The  intrinsic  capability  of  large-eddy  simulation  (LES)  for  predicting  sufficient  details 
of  unsteady  separating  flows  has  recently  been  explored  by  a  certain  number  of  researchers. 
Nagarajan  et  al.  (2006)  and  Ghias  et  al.  (2005)  performed  LES  of  dynamic  stall  over  a 
pitching  airfoil  and  tip- flow  of  a  rotor  in  hover  using  structured-grid  finite-difference  methods 
on  curvilinear  coordinates.  The  Reynolds  numbers  considered  in  these  simulations  were 
around  10°,  which  are  substantially  below  that  of  a  typical  helicopter  rotor  during  low- 
speed  maneuvers  (O(106)).  LES  at  higher  Reynolds  numbers  has  been  difficult  mainly 
due  to  the  numerical  instability  issue  and  high  computational  costs  associated  with  the 
spatial  and  temporal  resolution  requirements.  Simulations  and  experiments  performed  at 
Rec  <  5  x  105  indicate  that  stall  frequently  occurs  when  laminar  flow  separates  near  the 
leading  edge.  This  process  leads  to  large  steady-state  stall  hysteresis.  However,  turbulent 
separation  is  more  common  at  high  Reynolds  number  and  this  makes  the  unsteady  stall 
characteristics  quite  dissimilar  from  the  features  observed  at  low  Reynolds  number. 

Most  CFD  works  so  far  have  focused  on  the  validation  of  CFD  codes  and  qualitative 
features  of  dynamic  stall  rather  than  the  understanding  of  the  flow  physics  under  realistic 
flight  conditions.  Detailed  quantitative  characterization  of  velocity  and  pressure  are  needed 
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in  the  separating  flow  region  for  the  fundamental  understanding  of  the  unsteady  separation 
process.  In  this  work,  we  perform  wall-resolved  LES  of  unsteady  separation  over  a  pitching 
airfoil  at  realistic  Mach  and  Reynolds  numbers.  The  research  focus  is  on  a  quantitative 
understanding  of  the  unsteady  separation  process  rather  than  validation  and  qualitative 
characterization  of  the  flow  field.  For  this  purpose,  we  will  employ  an  unstructured-grid 
LES  technology,  which  maintains  numerical  stability  by  obeying  higher-order  conservation 
principles  —  i.e.,  kinetic-energy  conservation  in  the  inviscid  limit  in  addition  to  mass  and  mo¬ 
mentum  conservation.  The  unstructured  grid  topology  as  well  as  a  hybrid  imp  licit /explicit 
time-integration  method  provides  highly  enhanced  efficiency  in  treating  spatial  and  tempo¬ 
ral  resolution  requirements  in  the  dynamically  important  separated  flow  region. 

The  specific  objectives  of  the  proposed  work  are  (i)  to  develop  and  provide  guidelines 
for  an  accurate  LES  of  dynamic  stall  at  realistic  Reynolds  and  Mach  numbers  (Rec  = 
2  —  4  x  106,  Ma  =  0.2  —  0.4);  (ii)  to  gain  a  deeper  understanding  of  quantitative  aspects 
of  unsteady  separation  during  dynamic  stall,  such  as  turbulent  kinetic  energy  and  vorticity 
budgets,  velocity  and  pressure  fluctuations,  and  spatiotemporal  correlations  of  primitive 
variables;  (iii)  to  quantify  the  effects  of  the  important  parameters  on  dynamic  stall  such 
as  compressibility,  type  of  pitching  motion,  and  pitch  rate;  (iv)  to  provide  insights  for 
advancing  turbulence  models  by  analyzing  the  simulation  databases;  and  (v)  to  develop 
strategies  for  unsteady  separation  control. 
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Chapter  2 


Literature  survey 

Flow-field  visualization  has  been  a  dominant  method  to  characterize  or  categorize  the  dy¬ 
namic  stall  phenomenon.  For  example,  using  an  oil-smoke  visualization  technique,  Mc- 
Croskey  et  al.  (1976)  described  three  different  types  of  boundary-layer  separation:  (i) 
trailing-edge  stall;  (ii)  leading-edge  stall  following  a  progression  of  flow  reversal  from  the 
trailing  edge;  and  (iii)  leading-edge  stall  due  to  abrupt  bursting  of  a  leading-edge  separation 
bubble.  Carr  et  al.  (1991)  employed  a  real-time  point-diffraction  interferometry  technique 
to  obtain  interferograms  which  visualize  the  effect  of  unsteadiness  on  the  leading-edge  flow 
over  a  pitching  airfoil. 

Lorber  &  Carta  (1987)  was  probably  the  first  to  provide  an  experimental  database 
of  constant-pitch-rate  aerodynamic  information,  such  as  the  surface  pressure  and  overall 
loads,  at  realistic  combinations  of  Reynolds  and  Mach  numbers.  They  considered  an  airfoil 
oscillating  at  pitch  rates,  between  0.001  and  0.020,  Mach  numbers  between  0.2  and  0.4, 
and  Reynolds  numbers  between  2  —  4  x  106.  The  results  demonstrated  the  influence  of  the 
leading-edge  vorticity  on  the  unsteady  aerodynamic  response  during  and  after  stall.  The 
leading-edge  vortex  was  found  to  be  strengthened  by  increasing  the  pitch  rate  while  the 
vortex  strength  was  found  to  be  weakened  by  increasing  the  Mach  number  and  by  starting 
the  motion  close  to  the  steady-state  stall  angle.  In  the  experiment,  they  observed  a  periodic 
pressure  oscillation  after  stall  at  a  high  pitch  angle  and  moderate  Reynolds  number.  A  small 
supersonic  zone  near  the  leading  edge  at  Ma  =  0.4  was  found  to  reduce  significantly  the 
peak  suction  pressures  and  the  unsteady  increments  to  the  airloads.  Most  recently,  Sahoo 
et  al.  (2009)  conducted  particle  image  velocimetry  (PIV)  measurements  of  the  velocity  and 
pressure  during  the  initiation  of  dynamic  stall  processes  for  a  pitching  airfoil  at  realistic 
helicopter  flight  conditions.  Unlike  other  previous  experimental  work,  this  research  aimed 
at  gaining  a  quantitative  understanding  of  the  vorticity  and  turbulent  flow  characteristics 
in  the  separated  flow  region. 

Among  physical  processes  during  dynamic  stall,  the  onset  of  flow  separation  has  been 
studied  relatively  extensively.  Smith  (1988)  examined  several  important  issues  including  the 
instability  of  a  leading-edge  separation  bubble  and  finite-time  breakup  of  a  boundary  layer. 
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Later,  Currier  &  Fung  (1992)  performed  a  detailed  analysis  of  experimental  data  provided 
by  McCroskey  et  al.  (1982)  and  Lorber  &  Carta  (1987)  to  characterize  the  onset  of  stall 
for  airfoils  undergoing  oscillatory  and  ramp-type  motions  about  static  stall  angles.  They 
found  that,  although  the  dependence  of  stall  onset  on  the  Mach  number  is  the  same  for 
different  airfoils,  the  dependence  on  the  frequency  can  be  very  different  for  different  airfoils, 
or  for  airfoils  with  differently  shaped  leading  edges.  Currier  &  Fung’s  analysis  showed  that 
dynamic  stall  cases  considered  by  McCroskey  et  al.  (1982)  and  Lorber  &  Carta  (1987) 
belong  to  the  range  of  Reynolds  number  for  which  stall  onset  is  ascribed  to  bursting  of  a 
separation  bubble,  arising  from  failure  of  the  separated  boundary  layer  to  reattach.  It  was 
also  suggested  that  completion  of  the  separation  process  requires  a  finite  time  comparable 
to  other  flow  time  scales  and  that,  although  increased  unsteadiness  delays  the  occurrence  of 
the  onset  to  a  higher  angle,  it  actually  promotes  the  process  of  boundary-layer  separation. 

Much  research  has  also  been  conducted  to  identify  the  effects  of  key  parameters  such  as 
the  Mach  number,  Reynolds  number,  and  pitch  rate  on  the  dynamic  stall  phenomenon.  In 
wind-tunnel  experiments,  Jumper  et  al.  (1987)  reported  that  the  magnitude  of  the  maximum 
lift  coefficient  of  a  pitching  airfoil  is  highly  dependent  on  the  pitch  rate.  McCroskey  et  al. 
(1976)  investigated  the  influence  of  various  airfoil  profiles  and  leading  edge  geometries.  The 
effects  of  Mach  number  on  the  dynamic  stall  of  an  oscillating  airfoil  have  been  discussed  by 
Chandrashekhara  &  Carr  (1989). 

Although  these  and  other  previous  studies  have  provided  insights  into  qualitative  aspects 
of  the  breakdown  of  a  boundary  layer  and  onset  of  dynamic  stall,  and  some  effects  of  key 
operation  parameters,  quantitative  aspects  of  unsteady  separation  during  dynamic  stall, 
such  as  turbulent  kinetic  energy  and  vorticity  budgets,  velocity  and  pressure  fluctuations, 
and  spatiotemporal  correlations  of  primitive  variables,  are  yet  not  well  understood.  Such 
an  understanding  is  crucial  especially  for  advancing  theoretical  and  computational  methods 
and  models  for  predicting  dynamic  stall  and  for  developing  new  concepts  for  stall  control. 

Although  computational  fluid  dynamics  has  received  increasing  attention  as  a  means  to 
develop  an  understanding  of  detailed  spatiotemporal  characteristics  of  dynamic  stall  events, 
most  computational  studies  so  far  have  focused  on  the  validation  of  numerical  methods 
and  qualitative  features  of  dynamic  stall.  Visbal  (1988)  performed  RANS  simulations  of 
dynamic  stall  of  an  airfoil  which  was  pitched  at  a  constant  rate  from  zero  incidence  to 
a  high  angle  of  attack.  The  computed  dynamic  stall  events,  as  well  as  the  complicated 
effects  of  pitch  rate  and  axis  location  were  found  in  qualitative  agreement  with  experimental 
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observations.  He  also  investigated  compressibility  effects  on  dynamic  stall  and  found  that 
the  shock/boundary-layer  interaction  at  high  Mach  numbers  significantly  alters  the  dynamic 
stall  process  (Visbal  1990).  Later,  Visbal  (1991)  presented  a  computational  assessment  of 
various  techniques  for  dynamic  stall  control  using  RANS.  Choudhuri  et  al.  (1994)  conducted 
a  low  Reynolds  number  {Rec  =  104)  simulation  of  initial  stages  of  two-dimensional  unsteady 
leading-edge  boundary-layer  separation  of  laminar  subsonic  flow  over  a  pitching  NACA-0012 
airfoil  using  the  compressible  laminar  Navier-Stokes  equations. 

Large-eddy  simulation  (LES)  has  emerged  as  a  tool  for  studying  detailed  physics  of 
unsteady  separating  flows.  Ghias  et  al.  (2004,  2005)  performed  LES  of  dynamic  stall  over 
a  pitching  airfoil  and  tip-flow  of  a  rotor  in  hover  using  structured-grid  finite-difference 
methods  on  curvilinear  coordinates  at  relatively  low  Reynolds  numbers  ( Rec  =  104  and 
4.5  x  104).  They  could  successfully  compare  their  computational  results  with  the  experi¬ 
mental  data  of  Walker  et  al.  (1985).  They  could  also  utilize  the  LES  database  to  develop 
control  strategies  for  dynamic  stall.  Nagarajan  et  al.  (2006)  demonstrated  advantages  of 
LES  for  dynamic  stall  applications  by  conducting  comparative  RANS  and  LES  of  dynamic 
stall  over  a  pitching  NACA  0012  airfoil  at  a  transitional  Reynolds  number  ( Rec  =  1.3x  105). 
They  also  suggested  that  the  predictive  capability  of  RANS  not  only  for  the  flow-field  but 
also  for  the  radiated  noise  can  be  significantly  improved  with  refined  turbulence  models 
whose  construction  can  be  deduced  from  an  analysis  of  an  LES  database. 

The  Reynolds  numbers  considered  in  the  previous  LES  were  substantially  below  that  of 
a  typical  helicopter  rotor  during  low-speed  maneuvers  (O(106)).  This  was  mainly  because 
LES  at  higher  Reynolds  numbers  has  been  difficult  mainly  due  to  the  numerical  instability 
issue  and  high  computational  costs  associated  with  the  spatial  and  temporal  resolution 
requirements.  While  the  issue  of  high  computational  costs  may  be  relieved  by  alternative 
LES-related  methods  such  as  hybrid  LES-RANS  methods,  detached-eddy  simulation  (DES), 
and  kinetic-eddy  simulation  (KES)  (Fang  <fc  Menon  2006),  numerical  instability  remains  a 
major  obstacle  for  an  accurate  simulation  of  dynamic  stall  at  high  Reynolds  number. 
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Computational  methodology 

3.1  Flow  configuration 

The  flow  configuration  is  shown  in  figure  3.1  and  corresponds  to  the  experimental  setup 
of  Lorber  &  Carta  (1987),  which  was  developed  to  study  dynamic  stall  penetration  at 
constant  pitch  rate  and  realistic  combinations  of  Reynolds  and  Mach  numbers.  The  flow 
configuration  models  conditions  occurring  during  aircraft  post-stall  maneuvers  and  during 
helicopter  high  speed  forward  flight.  A  Sikorsky  SSC-A09  airfoil  with  a  chord  length  of 
43.9  cm  was  installed  in  the  UTRC  Wind  Tunnel.  The  surface  pressure  was  measured  using 
miniature  transducers,  and  the  locations  of  transition  and  separation  were  determined  using 
surface  hot  film  gages. 

There  are  four  basic  parameters  that  characterize  the  experimental  conditions:  Mach 
number,  Reynolds  number,  pitch  rate,  and  pitch  range.  While  the  majority  of  the  data  was 
taken  at  Mach  number  of  0.2,  some  data  at  Mach  numbers  of  0.3  and  0.4  were  also  recorded. 
These  Mach  numbers  cover  the  range  expected  both  for  aircraft  maneuvers  at  high  angle 
of  attack  and  helicopter  blades  on  the  retreating  side  of  the  rotor  at  high  advance  ratio. 
The  corresponding  Reynolds  numbers  based  on  the  airfoil  chord  were  2,  3,  and  4  x  106, 
respectively.  The  pitch  angle  was  controlled  by  hydraulic  actuators  attached  to  each  end  of 
the  airfoil.  Unsteady  oscillating  configurations  with  9  different  sinusoidal  and  36  different 
constant  pitch  rate  ramping  motions  were  considered  at  angles  of  attack  in  a  range  between 
—5°  and  30°.  Both  positive  and  negative  pitch  rate  ramps  were  studied  at  nondimensional 
pitch  rates,  A  =  dc/2C/00  =  0.001  —  0.02,  where  a,  c,  and  UQ 0  are  the  pitch  rate,  chord,  and 
freestream  velocity,  respectively.  For  the  sinusoidal  oscillations,  data  were  taken  at  reduced 
frequencies,  k  =  uc/2Uoo  =  0.025,  0.050,  and  0.100,  where  ui  is  the  angular  frequency. 

Pressure  time  histories,  aerodynamic  forces  and  moment,  surface  pressure  distributions, 
and  transition  and  separation  locations  are  available  as  a  function  of  pitch  rate,  type  of 
pitching  motion  (ramp  and  sinusoidal),  Mach  number,  and  Reynolds  number. 
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Figure  3.1.  Flow  configuration  for  LES  of  flow  over  a  pitching  airfoil. 

3.2  Numerical  methods 

We  employ  an  unstructured-grid  LES  solver  which  was  developed  at  the  Center  for  Tur¬ 
bulence  Research  (You  et  al.  2008a;  Shoeybi  et  al.  2007)  and  has  recently  been  further 
developed  by  the  PI  to  include  a  new  subgrid-scale  LES  model  (You  &  Moin  2007).  The 
numerical  method  is  based  on  unstructured-grid  finite-volume  discretization  of  the  Favre- 
filtered  compressible  Navier-Stokes  equations  with  subgrid-scale  stress  and  heat  flux  models. 
The  present  numerical  method  overcomes  two  major  difficulties  encountered  in  the  previous 
rotor  applications  using  structured-curvilinear-grid  LES  approaches  such  as  employed  by 
the  teams  at  the  Center  for  Turbulence  Research  (Nagarajan  et  al.  2006)  and  at  George 
Washington  University  (Ghias  et  al.  2005). 

Firstly,  it  is  known  that  non-  or  low-dissipative  finite-difference  schemes  on  curvilinear 
coordinates  do  not  strictly  conserve  kinetic  energy  (Nagarajan  et  al.  2004),  and  therefore,  are 
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more  prone  to  be  numerically  unstable  (see  figure  3.2(a)).  Due  to  the  numerical  instability, 
the  Reynolds  numbers  in  the  previous  LES  (Nagarajan  et  al.  2006;  Ghias  et  al.  2005)  were 
lower  by  an  order  of  magnitude  (O(105))  than  the  practical  Reynolds  number  in  rotor 
applications  (>  O(10b)).  The  present  Cartesian-coordinate-based  finite-volume  method 
maintains  the  numerical  stability  as  well  as  the  numerical  accuracy  at  high  Reynolds  number 
by  employing  an  unstructured-grid  spatial-discretization  algorithm  that  minimizes  the  non¬ 
conservation  of  kinetic  energy  (figure  3.2(6)). 

The  Favre-filtered  compressible  Navier-Stokes  equations  can  be  written  as 

dU  d  Fj  dGj  ,  . 
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U  is  the  vector  of  the  Favre-filtered  conserved  variables  and  Fj  and  Gj  are  the  flux 
vectors  in  the  j-direction.  p,  p,  m,  and  E  denote  density,  pressure,  velocity  component, 
and  energy,  respectively,  a ij  and  qj  are  the  filtered  stress  tensor  and  heat  flux,  respectively. 
t'-Js  and  qbgs  are  the  subgrid-scale  stress  tensor  and  heat  flux,  respectively. 

Finite- volume  discretization  of  the  governing  equation  (3.1)  leads  to 


dUk 

dt 


f  3= 1 


rij  =  0, 


(3.2) 


where  is  the  volume  measure  of  a  volume-element  f lkl  and  Uk  is  the  state  variable 
vector  at  grid  point  k.  F?  and  Gj  are  the  flux  vectors  at  the  element-boundary  faces  dlij,. 
rij  denotes  the  face-normal  unit  vector. 

In  the  present  study,  the  convective  and  diffusive  fluxes  are  obtained  using  the  Summation- 
By-Parts  operators  (e.g.,  skew-symmetric  averaging  for  the  convective  flux;  see  also  Strand 
(1994))  which  guarantee  non-growing  positive  norms  of  primary  variables,  thereby  main¬ 
taining  numerical  stability.  The  present  method  is  proven  to  be  particularly  well  suited  for 
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predicting  subtle  separation  effects  in  turbulent  boundary  layers  at  high  Reynolds  number 
where  boundary  layer  energetics  play  a  crucial  role  (You  et  al.  2008a). 

In  addition,  the  proposed  LES  will  employ  the  dynamic  global- coefficient  subgrid-scale 
(SGS)  model  which  has  recently  been  developed  by  the  PI  (You  Sz  Moin  2007).  In  the 
dynamic  Smagorinsky  model  (Germano  et  al.  1991;  Moin  et  al.  1991),  which  has  widely 
been  used  in  LES,  the  model  coefficient  is  dynamically  determined  as  a  function  of  space 
and  time  using  the  scale-invariance  concept  and  the  local- equilibrium  hypothesis  ( i.e an 
equilibrium  between  the  subgrid-scale  dissipation  and  the  viscous  dissipation  at  the  same 
physical  location).  Although  the  dynamic  model  coefficient  vanishes  where  the  flow  is 
laminar  or  fully  resolved,  it  can  cause  numerical  instability  since  its  value  often  becomes 
negative  and/or  highly  fluctuates  in  space  and  time. 

To  overcome  the  deficiency  of  the  dynamic  Smagorinsky  model,  the  PI  developed  a  dy¬ 
namic  procedure  for  determining  the  model  coefficient  utilizing  a  global  equilibrium  between 
the  subgrid-scale  dissipation  and  the  viscous  dissipation  (You  Sz  Moin  2007).  In  this  ap¬ 
proach,  the  model  coefficient  is  globally  constant  in  space  but  varies  in  time,  and  it  still 
guarantees  zero  eddy  viscosity  in  the  laminar-flow  regions.  The  model  does  not  require 
any  ad  hoc  numerical  stabilization  or  clipping  operation  which  is  usually  necessary  in  the 
local-equilibrium  based  dynamic  models. 

Secondly,  it  has  been  difficult  to  selectively  resolve  dynamically  important  separated 
flow  regions  using  H-  and  O-type  curvilinear  mesh  topologies,  which  are  most  commonly 
used  in  rotor  applications.  Previous  experience  with  wall-resolved  LES  of  turbulent  flows 
indicates  that  streamwise  and  spanwise  spacings  of  about  50-100  and  30-50  wall-units, 
respectively,  are  required  in  the  separated  flow  region  and  wake.  Although  local-refinement 
or  overset  type  grids  may  be  used,  it  is  known  that  numerical  errors  associated  with  the 
spatial  interpolation  adversely  affect  the  simulation  results.  In  the  present  method,  the 
issue  is  overcome  by  an  unstructured-grid  topology,  which  provides  higher  flexibility  and 
efficiency  in  distributing  mesh  resolution  (figure  3.3(a)). 

Furthermore,  the  high  CFL  number  restriction  in  the  local  dense  mesh  region  is  allevi¬ 
ated  with  the  use  of  a  hybrid  imp  licit /exp  lict  time-advancement  scheme.  The  discretized 
governing  equation  (3.1)  is  recast  as 

dU 

=  H(U),  (3.3) 
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where  U  =  {Uf,  U{  ,  Uj  ,  •  •  •  ,  U'^r)T  is  the  solution  vector  and  N  is  the  number  of  grid  points. 
Similarly,  the  right-hand  side  H  =  [HJ .  Hf,  H±  ,  •  •  •  ,  Hjj)T  is  the  flux  vector  defined  as 


Hk 


E  t{Fi 

f  3  =  1 


G{ 


f 

4 


(3.4) 


The  right-hand  side  of  equation  (3.3)  is  decomposed  as  H  =  He  +  Hl,  where  He  and  Hl 
are  the  explicit  and  implicit  parts  of  H .  The  stiffness  of  H  is  estimated  from  eigenvalues  of 
the  flux  Jacobian  matrix  J  as  follows 


J=w-W 


in  -  BJIi 

•J  ‘ATT  ' 

dUj 


z  —  1,  2, 3,  •  •  •  ,  N,  j  —  1,  2,  3,  •  •  •  ,  N. 


Since  the  computation  of  eigenvalues  is  very  costly  for  large-scale  simulations,  we  utilize 
the  Gerschgorin  theorem  which  gives  a  radius  of  a  disc  including  all  the  eigenvalues  of  the 
Jacobian  matrix,  thereby  providing  maximum  allowable  time  step  size  for  stable  zone  of 
an  explicit  scheme.  The  decomposed  equation  is  advanced  in  time  using  a  semi-implicit 
Runge-Kutta  method  proposed  by  Le  &  Moin  (1991). 

This  feature  allows  to  adaptively  divide  the  mesh  into  explicit  and  implicit  zones  (fig¬ 
ure  3.3(6))  and  leads  to  significantly  reduced  memory  requirements  while  maintaining  the 
advantage  of  an  implicit  integration  scheme.  This  method  is  especially  advantageous  for 
long-time  integration  of  low-pitch-rate  unsteady  separating  flow. 
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Figure  3.2.  Contours  of  the  instantaneous  streamwise  velocity  around  the  leading-edge  of  an 
airfoil  predicted  by  (a)  LES  based  on  a  second-order  centered  finite-difference  method  on  curvilinear 
coordinates  and  by  (b)  LES  based  on  a  second-order  centered  finite-volume  method  on  Cartesian 
coordinates  (the  present  scheme). 
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Figure  3.3.  Schematic  illustration  of  (a)  a  grid  resolution  topology  and  (b)  separation  of  implicit 
and  explicit  time-integration  zones. 
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Chapter  4 


Flow  over  a  pitching  airfoil 

4.1  Effects  of  mesh  resolution 

A  series  of  computational  grids  were  employed  to  investigate  the  effects  of  mesh  resolution 
on  the  prediction  of  flow  over  a  steady  and  pitching  airfoils.  As  discussed  in  section  3.2,  the 
present  method  utilizes  advantages  of  using  unstructured  grid  and  adaptive  implicit-explicit 
time-integration  scheme.  Each  grid  consists  of  multiple  zones  with  different  grid  resolution. 
Grid  lines  are  clustered  around  the  suction  surface  of  the  airfoil  and  in  the  wake  region  so 
that  the  separated  shear  layer,  recirculating  flow,  and  turbulent  wake  are  well  resolved. 

Most  of  the  computational  domain  is  discretized  using  hexahedral-shape  elements.  Es¬ 
pecially  near  the  airfoil  surface,  an  elliptic-type  mesh  generator  is  employed  to  align  grid 
lines  to  be  parallel  and  orthogonal  to  the  airfoil  surface  in  the  streamwise  and  wall-normal 
directions,  respectively.  As  illustrated  in  figure  4.1,  there  are  numbers  of  interfaces  where 
mesh  resolution  is  transitioned  from  fine  to  coarse.  These  interfaces  are  patched  with  prism- 
shape  elements  to  maintain  the  quality  of  arrangement  of  primitive  variables. 

A  mesh  with  about  6  million  cells  was  designed  first  and  employed  for  a  coarse-resolution 
LES  of  flow  over  a  steady  airfoil  at  a  fixed  angle  of  attack  in  order  to  assess  the  resolution 
requirement.  Subsequently,  a  24  million  cell  mesh  is  designed. 

Figure  4.2  shows  pressure  distributions  on  the  surface  of  the  Sikorsky  SSC-A09  airfoil 
at  a  fixed  angle  of  attack  of  14  degrees.  The  Reynolds  number  and  Mach  number  are 
2  x  106  and  0.2,  respectively.  The  pressure  distribution  predicted  by  the  present  LES  on  a 
24  million-cell  mesh  is  found  to  agree  well  with  the  experimental  measurement  in  Lorber  <fc 
Carta  (1987).  The  near  wall  resolution  in  wall  units  on  the  24  million-cell  mesh  is  found  to 
be  A.x+  =  50  -  450,  A y+  =  1-2,  and  Az+  =  60  -  130. 

4.2  Results  and  discussion 

The  24  million  cell  mesh,  which  was  demonstrated  to  be  reasonably  capable  of  predicting 
flow  over  the  Sikorsky  SSC-A09  airfoil  at  the  operating  Reynolds  and  Mach  numbers,  is  also 
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employed  for  large-eddy  simulations  of  flow  during  a  pitching  motion  of  the  airfoil.  To  avoid 
mesh  rotation  or  re-meshing  during  the  pitching  motion,  the  Favre-filtered  compressible 
Navier-Stokes  equations  are  recast  into  forms  in  a  non-inertial  reference  frame.  Therefore, 
large-eddy  simulations  are  performed  for  variables  in  a  non-intertial  reference  frame  along 
with  the  effects  of  the  Coriolis,  centrifugal,  and  rotational  acceleration  forces. 

Simulations  are  ongoing  for  two  different  cases.  The  first  case  corresponding  to  an 
experimental  case  conducted  at  Mach  number  of  0.2  and  a  sinusoidal  pitching  motion  of 
which  angle  of  attack  varies  in  time  as  follows: 

a(t)  =  20°  —  10°cos(cuf), 

where  a  is  the  angle  of  attack  as  defined  in  figure  3.1  and  u  is  the  reduced  frequency  and 
is  0.10  x  2Uoo/c.  The  second  case  also  corresponds  to  an  experimental  case  conducted  at 
Mach  number  of  0.3  and  a  sinusoidal  pitching  motion  as  follows: 

a(t)  =  12°  —  8°  cos  (cot), 

where  the  reduced  frequency  is  identical  to  the  first  case. 

These  two  simulations  are  ongoing  while  the  present  simulations  confirm  the  stability  and 
effectiveness  of  the  presented  numerical  schemes  for  dynamic  stall  simulations  at  realistic 
operating  Reynolds  and  Mach  numbers. 

Large-eddy  simulations  are  being  conducted  with  the  acoustic  Courant-Friedrichs-Lewy 
(CFL)  numbers  of  25  and  20  for  Mach  number  of  0.2  and  0.3  cases,  respectively.  The 
acoustic  CFL  numbers  correspond  to  the  time-step  sizes  of  0.27  x  10-3  and  0.23  x  10~3 
normalized  by  the  airfoil  chord  length  and  the  speed  of  sound  for  Mach  0.2  and  0.3  cases, 
respectively.  Using  512  cores  of  SGI  Altix  ICE  8200LX  computer,  about  12  -  16  days  are 
required  for  large-eddy  simulation  over  a  pitching  period  when  a  24  million-cell  mesh  is 
employed. 

Figures  4.3  and  4.4  show  gross  features  of  flow  over  pitching  airfoils  at  two  different  con¬ 
ditions.  The  characteristics  of  flow  separation  and  reattachment  processes  are  qualitatively 
congruent  with  experimental  observation  by  Lorber  &;  Carta  (1987). 

In  figure  4.5,  the  lift  coefficient  and  the  pressure  drag  coefficient  over  the  pitching  airfoil 
are  compared  against  the  experimental  data  of  Lorber  &  Carta  (1987)  as  a  function  of  pitch 
angle  during  a  pitching  cycle  at  Mach  number  of  0.2.  Overall,  the  lift  coefficient  is  in  good 
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agreement  with  the  experimental  data.  However,  in  the  early  stage  of  downward  pitching 
motion  (30°  — >•  23° ) ,  a  noticeable  deviation  of  the  pressure  drag  is  predicted  by  the  present 
LES  when  compared  to  the  experimental  data. 

The  pitching  moment  coefficient  over  a  pitching  cycle  is  plotted  in  figure  4.6  (top)  and 
is  found  to  well  agree  with  the  experimental  measurement.  The  variation  of  the  pressure 
coefficient  near  the  leading  edge  of  the  pitching  airfoil  (x/c  =  0.005)  is  also  found  to  be  in 
good  agreement  with  the  experimental  data  as  shown  in  figure  4.6  (bottom). 

4.3  Summary 

To  gain  a  quantitative  understanding  of  the  unsteady  separation  process  over  a  pitching 
airfoil,  large-eddy  simulations  (LES)  of  turbulent  flow  over  a  pitching  airfoil  at  realistic 
Reynolds  and  Mach  numbers  have  been  conducted.  A  novel  combination  of  discretization 
schemes  and  time-integration  schemes  was  employed  to  achieve  numerical  stability  at  high 
Reynolds  number  simulations  through  higher-order  conservation  and  to  be  equipped  with  a 
highly  efficient  way  to  treat  time-step  size  restriction  in  the  separated  flow  region  locally  re¬ 
fined  with  dense  mesh.  The  present  simulations  confirm  the  stability  and  effectiveness  of  the 
presented  numerical  schemes  for  dynamic  stall  simulations  at  realistic  operating  Reynolds 
and  Mach  numbers  and  show  the  characteristics  of  flow  separation  and  reattachnrent  pro¬ 
cesses  which  are  qualitatively  congruent  with  experimental  observation. 
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Figure  4.1.  Computational  grid  topology  for  LES  of  flow  over  a  pitching  airfoil. 
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x/c 

Figure  4.2.  Pressure  distribution  on  the  airfoil  surface  at  the  angle  of  attack  14  degrees.  Solid 
line,  the  present  LES  solution  on  a  24  million-cell  mesh;  symbol,  experimental  data. 
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20° 
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Figure  4.3.  Contours  of  the  spanwise  vorticity  predicted  by  the  present  LES  over  a  half  period  of 
pitching  motion  at  Mach  number  of  0.2. 
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Figure  4.4.  Contours  of  the  spanwise  vorticity  predicted  by  the  present  LES  over  a  period  of 
pitching  motion  at  Mach  number  of  0.3. 
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Figure  4.5.  Lift  coefficient  and  pressure  drag  coefficient  as  a  function  of  pitch  angle  at  Mach 
number  of  0.2. 
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time/period 


Figure  4.6.  Pitching  moment  coefficient  and  pressure  coefficient  at  x/ c  =  0.005  as  a  function  of 
pitch  angle  at  Mach  number  of  0.2. 
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Chapter  5 


Direct  numerical  simulation  of  unsteady 
separation  of  turbulent  boundary  layers 

5.1  Introduction 

The  term,  turbulent  flow  separation,  describes  a  process  of  detachment  of  turbulent  boundary- 
layer  flow  from  the  wall.  In  contrast  to  a  steady  separating  turbulent  boundary  layer,  where 
the  detachment  location  is  statistically  fixed,  in  an  unsteady  separating  turbulent  bound¬ 
ary  layer,  the  flow  detachment  location  varies  in  time  due  to  an  organized  time-dependent 
flow  condition  and/or  a  time- varying  wall  motion  (Simpson,  1989;  McCroskey,  1982).  Un¬ 
derstanding  unsteady  separation  of  turbulent  boundary  layers  has  considerable  practical 
significance  because  unsteady  separation  often  characterizes  the  performance  and  efficiency 
of  many  aerodynamic  applications  such  as  helicopter  rotor  blades,  wind  turbine  blades, 
pitching  and  flapping  airfoils  and  wings,  and  rotating  turbomachinery  blades. 

For  instance,  helicopter  rotor  blades  experience  a  pitching  motion  which  leads  to  un¬ 
steady  flow  separation  followed  by  load  and  pitching-moment  overshoots.  The  unsteady 
flow  separation  can  in  turn  lead  to  unacceptably  large  vibratory  loads  and  acoustic  noise, 
and  limit  forward  flight  speeds  and  maneuverability.  Steady  separation  of  a  turbulent 
boundary  layer  over  an  airfoil  is  often  converted  to  unsteady  separation  when  the  separa¬ 
tion  is  controlled  by  oscillatory  blowing-suction  actuations  such  as  synthetic  jets  (Gilarranz 
et  al. ,  2005).  The  separation  point  over  a  composite-material  propeller  blade  becomes  time- 
dependent  as  the  blade  is  deformed  during  the  operation. 

In  spite  of  past  experimental  and  numerical  studies  on  separating  flows  in  both  modeled 
and  realistic  configurations,  quantitative  aspects  of  unsteady  separating  turbulent  bound¬ 
ary  layers,  such  as  Reynolds  stress  and  vorticity  budgets,  velocity  and  pressure  fluctuations, 
and  spatiotemporal  correlations  of  primitive  variables,  are  not  well  understood.  Agarwal 
&  Simpson  measured  phase-averaged  velocity  profiles  in  an  unsteady  separating  turbulent 
boundary  layer  over  a  flat  plate  where  the  unsteady  adverse  pressure  gradient  was  imposed 
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by  time- varying  upper- wall  suction  (Agarwal  &  Simpson,  1990).  They  suggested  that  quan¬ 
titative  measurements  of  Reynolds  stresses  and  turbulence  structures  near  the  wall  in  the 
separated  flow  region  are  necessary  to  improve  our  understanding  of  the  unsteady  separa¬ 
tion  process.  However,  most  experimental  studies  have  concentrated  on  measurements  of 
aerodynamic  forces  such  as  the  surface  pressure  and  overall  loads,  or  on  the  flow- field  visual¬ 
ization.  Quantitative  measurements  of  the  unsteady  separated  flow  field  have  been  difficult 
using  experimental  techniques,  and  therefore,  have  rarely  been  reported  in  the  literature. 

Computational  works,  especially  at  practical  Reynolds  numbers,  have  mostly  been  per¬ 
formed  using  the  Reynolds-averaged  Navier-Stokes  (RANS)  equations  or  its  unsteady  coun¬ 
terpart  (URANS)  (e.g.,  Visbal  (1988);  Spentzos  et  al.  (2005)).  However,  it  is  known  to  be 
very  challenging  for  (U)RANS  to  accurately  predict  highly  unsteady  flow  involving  incipient 
flow  separation,  formation  and  evolution  of  stall  vortices,  and  reattachment.  So  far,  none 
of  the  available  turbulence  models  have  shown  a  satisfactory  predictive  performance  for  all 
of  these  flow  phenomena.  A  comprehensive  database,  especially  Reynolds  stress  equation 
budgets,  would  be  highly  valuable  for  progress  in  this  area.  In  a  review  article,  Carr  & 
McCroskey  (1992);  McCroskey  (1982)  concluded  that  under  this  condition  (unsteady  sep¬ 
aration),  the  use  of  a  turbulence  model  based  on  equilibrium  attached  boundary  layers  in 
steady  flow  {e.g.,  eddy  viscosity,  Baldwin-Lomax)  is  open  to  question.  The  task  of  predicting 
separation  by  definition  deals  with  boundary  layers  that  have  experienced  strong  pressure 
gradients,  often  both  positive  and  negative;  the  flow  approaching  unsteady  separation  con¬ 
tains  high  levels  of  vorticity  induced  by  these  pressure  gradients,  and  is  strongly  unsteady. 
Recent  study  has  shown  that  modification  of  the  turbulence  model  can  completely  change 
the  resultant  flow;  at  the  same  time,  very  little  has  been  experimentally  documented  about 
the  character  of  turbulence  under  these  conditions. 

The  unsteady  separation  process  is  known  to  be  highly  dependent  on  the  unsteady 
characteristics  and  strength  of  the  adverse-pressure  gradient  while  the  dependence  on  the 
Reynolds  number  is  virtually  unknown  (McCroskey,  1982).  Although,  in  many  cases,  the 
surface  curvature  also  influences  the  separation  of  turbulent  boundary  layers,  it  would  be 
instructive  to  isolate  the  effects  of  the  unsteady  adverse  pressure  gradient  by  considering 
a  flat  plate  surface  with  imposed  adverse  pressure  gradients.  Coleman  &  Spalart  (1993) 
and  Na  &;  Moin  (1998)  performed  direct  numerical  simulations  (DNS)  of  steady  separated 
turbulent  boundary  layers  on  a  flat  plate  with  an  imposed  steady  adverse  pressure  gradient. 
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They  were  able  to  provide  insight  into  the  characteristics  of  steady  separated  flows  and  flow- 
field  databases  containing  complete  budgets  of  Reynolds  stresses.  However,  for  unsteady 
separating  turbulent  boundary  layers,  such  detailed  characterization  of  flow  physics  and 
comprehensive  flow-field  databases  are  not  available  so  far. 

In  this  study,  we  will  employ  a  computational  setup  consisting  of  a  turbulent  boundary 
layer  over  a  flat  plate  with  a  deliberate  unsteady  adverse  pressure  gradient  imposed  by 
time- varying  blowing-suction  on  the  upper  boundary.  The  setup  will  allow  us  to  focus  on  a 
quantitative  physical  understanding  of  the  unsteady  separation  process  including  boundary- 
layer  detachment  and  reattachment,  and  production  and  dissipation  of  turbulent  kinetic 
energy  and  vorticity.  An  understanding  of  the  characteristics  of  pressure  fluctuations  in 
unsteady  separating  turbulent  boundary  layers  is  also  anticipated.  Since  the  Reynolds 
shear  stress  and  its  gradient  are  largest  away  from  the  wall  in  separated  flows,  it  has 
been  suggested  that  the  largest  pressure  fluctuations  are  not  at  the  wall  but  away  from  it 
(Gilarranz  et  ai,  2005).  Identifying  the  structure  of  pressure  fluctuations  away  from  the 
wall  is  critical  to  the  prediction  and  understanding  of  turbulent  boundary  layer  acoustics. 
DNS  is  ideal  for  the  analysis  of  pressure  fluctuations  because  static  pressure  within  the  flow 
cannot  be  measured  experimentally. 

5.2  Computational  methodology 

5.2.1  Flow  configuration 

The  present  DNS  configuration  is  motivated  by  the  experimental  configuration  of  Agarwal 
&  Simpson  (1990)  and  is  schematically  shown  in  figure  5.1.  In  the  present  DNS,  the  in¬ 
flow  is  a  fully-developed  zero-pressure-gradient  turbulent  boundary  layer  and  is  generated 
using  a  recycling  technique  developed  by  Lund  et  al.  (1998).  No-slip  boundary  conditions 
are  used  along  the  lower  boundary  of  the  computational  domain  while  periodic  boundary 
conditions  are  applied  in  the  spanwise  direction.  At  the  exit,  convective  outflow  conditions 
are  employed  to  allow  turbulence  structures  to  smoothly  leave  the  computational  domain. 

A  novel  feature  of  the  present  computational  setup  is  that  deliberate  unsteady  adverse 
pressure  gradients  (e.g.,  replicating  the  oscillatory  adverse  pressure  gradient  in  Agarwal  & 
Simpson  (1990)’s  experiment)  can  be  imposed  by  Vtop(x,t)  which  represents  time-varying 
blowing  and  suction  at  the  upper  boundary.  One  way  of  obtaining  approximate  but  rea¬ 
sonable  VtoP(x,t )  is  to  use  a  potential  flow  solution  (from  the  panel  method,  for  example). 
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Figure  5.1.  Schematic  illustration  of  computational  modeling  of  unsteady  separating  turbulent 
boundary  layers. 


Several  iterations  of  viscous-inviscid  calculations  can  be  performed  to  improve  the  approx¬ 
imation.  At  the  upper  boundary,  the  streamwise  velocity  Utop(x,t)  is  adjusted  from  a 
zero-vorticity  condition  prescribed  by  the  equations  (5. 1,5. 2, 5. 3): 


v(x,Ly,z,t)  =  VtoP(x ), 


(5.1) 


du  dVt0p(x') 

dy  x,Ly,z,t  dx 


dw 

dy 


X,Ly,Z,t 


=  0, 


(5.2) 

(5.3) 


where  Ly  is  the  height  of  the  domain.  This  boundary  condition  assures  zero  vorticity  in  the 
spanwise  and  streamwise  directions. 

The  distinct  characteristics  of  unsteady  separating  turbulent  boundary  layers  are  re¬ 
vealed  by  a  systematic  comparison  with  steady  attached/separated  turbulent  boundary 
layers.  For  this  purpose,  four  different  flow  configurations  are  simulated:  (i)  turbulent 
boundary  layer  flow  with  a  zero-pressure  gradient  (ZPG);  (ii)  turbulent  boundary  layer 
flow  with  an  adverse-pressure  gradient  (APG);  (iii)  steady  separated  turbulent  boundary 
layer  flow  (SBL);  and  (iv)  unsteady  separating  turbulent  boundary  layer  flow  (USBL). 

For  all  configurations  grid  spacings  in  the  streamwise  and  spanwise  directions  are  uni¬ 
form.  In  the  wall-normal  direction,  the  grid  is  stretched  based  on  a  one-parameter  hyper¬ 
bolic  tangent  function  described  by  Na  &  Moin  (1998).  All  simulations  are  performed  in  a 
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configuration  similar  to  that  employed  by  Na  &  Moin  (1998).  The  extent  on  the  domain  in 
the  streamwise  direction  is  equivalent  to  350d*n,  where  S*n  is  the  displacement  thickness  at 
the  inlet  of  the  computational  domain.  The  inlet  location  is  defined  as  where  x  =  0.0,  and 
is  marked  by  the  dotted  line  in  figure  5.2.  The  spanwise  and  wall- normal  domain  sizes  are 
50  and  64<5*n,  respectively,  as  shown  schematically  in  figure  5.2.  The  Reynolds  number  is 
fixed  to  300,  based  on  the  momentum  thickness  and  freestream  velocity  at  the  inlet  in  all 
four  flow  configurations. 


Figure  5.2.  Computational  domain  sizes  in  terms  of  the  inlet-displacement  thickness  and  corre¬ 
sponding  number  of  grid  points. 

The  Vt0p(x)  prescribed  for  the  APG  and  SBL  cases  are  shown  in  figures  3  (a)  and  (b), 
respectively,  and  correspond  to  the  boundary  conditions  developed  by  Na  &  Moin  (1998). 
The  inflow  zero-pressure  gradient  turbulent  boundary  layer  is  generated  by  a  recycling 
technique  developed  by  Lund  et  al.  (1998)  in  the  upstream  region  which  extends  the  domain 
by  44 5*n,  (L'x  in  figure  5.2).  In  the  USBL  case  a  time-varying  Vtop(x,t),  which  is  similar  in 
shape  to  that  of  SBL  but  oscillates  in  amplitude,  is  employed  as  illustrated  in  figure  5.1. 

5.2.2  Numerical  methods 

The  incompressible  Navier-Stokes  equations  are  discretized  in  space  using  a  second-order 
central-difference  scheme  on  a  staggered  grid.  The  discretized  governing  equations  are  in¬ 
tegrated  in  time  using  a  semi-implicit  scheme.  A  low-storage  third-order  Runge-Kutta 
scheme  is  employed  for  treating  convective  terms  explicitly  and  a  second-order  Crank- 
Nicolson  scheme  is  employed  for  treating  viscous  terms  implicitly.  The  hybrid  Runge- 
Kutta/Crank-Nicolson  scheme  is  combined  with  a  modified  fractional-step  procedure,  where 
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Figure  5.3.  Vertical  velocity  profiles  imposed  on  the  top  boundary  of  the  computational  domain. 
Velocity  magnitudes  are  normalized  with  the  freestream  velocity. 
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the  divergence-free  velocity  field  is  obtained  by  solving  the  pressure  Poisson  equation  only 
at  the  last  substep.  The  Poisson  equation  is  solved  using  multi-grid  with  Fast  Fourier 
Transform  in  the  spanwise  direction.  The  numerical  algorithm  and  solution  methods  are 
described  in  detail  in  You  et  al.  (2007).  The  code  has  been  modified  for  simulations  of 
unsteady  separating  turbulent  boundary  layers.  The  present  solver  has  been  extensively 
validated  in  several  cases  of  internal  and  external  turbulent  flows  (e.g.,  You  et  al.  (2007)). 


5.3  Results  and  discussion 

5.3.1  Steady  boundary  layers 

The  present  DNS  results  in  the  steady  attached  (ZPG  and  APG)  and  separated  turbulent 
boundary  layer  (SBL)  configurations  are  extensively  validated  against  other  available  DNS 
data.  Data  collected  from  the  ZPG  configuration  agrees  very  well  with  Na  &  Moin’s  results 
(Na  &  Moin,  1998).  Some  results  of  ZPG  are  shown  in  figure  5.4(b)  and  compared  against 
results  by  Na  &  Moin  (1998).  Both  the  mean  velocity  and  root-mean-squared  fluctuating 
velocity  profiles  are  recorded  at  92 8*n  downstream  from  inlet. 

Data  collected  from  APG  are  shown  in  figures  5.5  and  5.6.  Velocity  profiles  in  wall- 
coordinates  and  root-mean-squared  profiles  are  shown  for  three  downstream  distances.  In 
the  upstream  location,  the  mean  velocity  and  rms  velocity  fluctuation  profiles  are  found  to  be 
nearly  identical  to  those  in  the  ZPG  case.  The  change  in  the  shape  of  the  wall-coordinate 
velocity,  shown  in  figure  5.5,  is  a  result  of  the  imposed  adverse  pressure  gradient.  The 
increase  of  the  magnitude  near  the  top  boundary  is  a  result  of  mass  being  inserted  into  the 
domain  in  the  suction  region  of  V)op(x).  The  effect  of  an  adverse  pressure  gradient  on  the 
distributions  of  the  root-mean-squared  profiles  is  shown  in  figure  5.6  in  three  streamwise 
locations.  With  the  increasing  streamwise  distance,  the  peak  of  the  streamwise  velocity 
fluctuation  profile  gradually  decreases  in  magnitude  while  the  peak  of  the  cross-stream 
velocity  fluctuation  profiles  increases  slightly,  illustrating  a  slight  increase  in  isotropy. 

Figure  5.7  shows  the  mean  velocity  profiles  in  the  SBL  case  in  several  upstream  and 
downstream  locations  from  the  separation  bubble.  The  time-averaged  location  of  the  sep¬ 
aration  bubble  spans  approximately  from  150  to  260  inlet-displacement  thicknesses  (S*n) 
downstream  of  the  inlet.  The  increase  in  the  streamwise  velocity  is  due  to  the  decrease  in 
wall  shear  stress  as  the  flow  approaches  the  separation  bubble  (figure  5.7(a)).  Unlike  in 
the  APG  configuration,  the  blowing  region  precedes  the  suction.  As  a  result,  the  mass  flux 
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(b)  RMS  velocity  fluctuations  normalized  with  the  freestream 
velocity 

Figure  5.4.  DNS  results  in  the  ZPG  case.  Na  &  Moin’s  results  shown  with  dotted  lines. 


34 


CHAPTER  5.  DNS  OF  UNSTEADY  SEPARATION 


Figure  5.5.  Mean  streamwise  velocity  profiles  in  the  APG  case  at  three  streamwise  locations  in 
wall  units. 

through  the  domain  is  minimized  at  the  inflection  point  of  Vtop{x),  approximately  220 5*n 
downstream  of  the  inlet.  Thus,  unlike  APG  results  (figure  5.5),  where  the  mean  velocity 
near  the  upper  boundary  increases  through  the  adverse  pressure  gradient,  it  actually  de¬ 
creases  in  the  SBL  configuration  (figure  5.7(b)).  Downstream  of  the  separation  bubble, 
the  shear  stress  begins  to  recover  its  original  magnitude.  This  recovery  does  not  require  a 
considerable  distance. 
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Profiles  of  the  root-mean-squared  velocity  fluctuations  in  the  SBL  case  are  shown  in 
figure  5.8.  Profiles  in  figures  5.8(a)-(c)  are  computed  inside  the  separation  bubble,  and 
profiles  in  figure  5.8(d)  are  computed  well  after  reattachment.  Figure  5.8(b)  shows  that 
the  maximum  turbulent  kinetic  energy  is  found  away  from  the  wall  where  the  separation 
bubble  height  is  maximized.  This  is  because  the  energy-containing  turbulent  structures  are 
lifted  around  the  separation  bubble.  Downstream  the  near-wall  peak  of  turbulent  kinetic 
energy  reappears,  however  the  overall  profiles  are  significantly  altered  from  those  in  zero- 
pressure-gradient  flow.  Figure  5.9  illustrates  the  skin-friction  coefficient  throughout  the 
computational  domain.  The  change  in  the  average  spanwise  spacing  of  the  streaks  from 
upstream  to  downstream  of  the  separation  bubble  is  clearly  identified. 

5.3.2  Unsteady  separating  turbulent  boundary  layer 
Flow  configuration 

As  mentioned  in  section  5.2.1  the  time  varying  adverse  pressure  gradient  is  imposed  by 
Vtop{x,t)  in  the  USBL  case.  The  Vtop(x )  used  in  the  SBL  case,  is  multiplied  by  an  ampli¬ 
fication  factor  which  is  a  function  of  time  to  create  Vtop(x,t)-  The  Vtop(x,t)  used  for  the 
present  USBL  case  employs  a  compound  trigonometric  function  which  produces  a  synco¬ 
pated  amplification  factor  ranging  between  1.0  and  0.6: 

VtoP(x,t)  =  Vtop(x)  *  [cos (cut  —  sin (ut))  +  4.0J/5.0  (5.4) 

The  period,  — ,  is  equivalent  to  1.5  flow  through  times  based  on  the  freestream  velocity  and 
the  potential  separation  bubble  length,  i.e.  the  maximum  distance  between  detachment 
and  reattachnrent  locations.  One  cycle  of  the  amplification  factor  as  a  function  of  the  phase 
number  is  shown  in  figure  5.10. 

Unsteady  flow-reversal 

The  time-dependent  boundary  condition  results  in  a  separation  bubble  which  intermittently 
stops  circulating,  and  subsequently  becoming  what  will  be  referred  to  as  a  hump.  Figure  5.11 
shows  the  spanwise-averaged  fraction  of  time  (yu)  that  the  flow  moves  downstream  for  a 
series  of  phases.  The  stalling  of  the  circulation,  and  full  attachment,  is  a  direct  result  of 
the  changes  in  streamwise  pressure  gradient.  It  should  be  noted  that  the  pressure  field  at 
a  given  time,  r,  is  not  only  a  function  of  the  magnitude  of  Vtop(x,T),  but  also  depends  on 
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its  time  rate,  The  phase  of  the  pressure  field,  and  circulation  response,  lead  the 

phase  of  the  amplification  factor. 

The  skin-friction  coefficient,  Cf,  depicts  the  region  in  which  flow  is  reversed.  Figure  5.12 

illustrates  the  streamwise  distribution  of  Cf  and  pressure  coefficient,  Cp,  evaluated  at  the 

wall.  Time-averaged  Vtop,max{x,t)  closely  resembles  that  in  the  SBL  case.  However,  the 

phase-averaged  pressure  distribution  in  each  phase  deviates  significantly  from  the  time- 

averaged  pressure  distribution.  The  pressure  held  is  shifted  such  that  a  point  near  the  top 

boundary  in  the  far  upstream  remains  zero  throughout  the  simulation.  The  variation  in 

wall  Cp,  |||  ,,,  leads  the  changes  in  Cf  by  approximately  one  phase.  In  the  SBL  case, 

the  plateaus  in  both  curves  correspond  closely  in  the  streamwise  location.  In  this  USBL 

dC 

case,  Cp  varies  considerably,  while  Cf  exhibits  little  shifting  in  location,  shows  little 
variation  in  the  boundary  layer  recovery  region,  x/5*n  >  300. 

Mean  velocity  field 

As  shown  in  figure  5.13,  the  averaged  skin-friction  and  pressure  coefficients  in  the  USBL 
case  are  found  to  be  nearly  identical  to  those  in  the  SBL  case.  This  suggests  that  strong 
similarities  would  exist  in  other  mean  flow  statistics  between  the  USBL  and  SBL  cases. 

Upstream  of  the  bubble  the  mean  streamwise  velocity  profiles  agree  well  each  other 
(figure  5.14(a)).  Inside  the  separation  bubble  in  the  SBL  case,  much  stronger  backflow 
is  observed  (figure  5.14(b)).  This  would  be  expected  because  during  the  USBL  case  the 
circulation  is  completely  seized  for  at  least  ten  percent  of  the  cycle  period.  As  a  result  of 
the  decrease  in  average  backflow  the  USBL  case  recovers  more  rapidly  (figure  5.14(c)).  And, 
downstream  of  the  bubble,  the  similarity  reappears  (figure  5.14(d)),  where  both  exhibit  a 
depression  below  the  log-law  region,  as  observed  by  Na  &;  Moin  (1998). 

Figures  5.15(a)-(d)  illustrate  differences  in  turbulence  intensities  between  the  USBL  and 
SBL  cases.  While  both  SBL  and  USBL  cases  agree  in  the  upstream  location  (figure  5.15(a)), 
the  maxima  in  the  SBL  are  shifted  more  away  from  the  wall  than  those  in  the  mean  USBL 
case  (figures  5.15(b)  and  (c)).  In  following  sections,  the  changes  in  the  separation  bubble 
height  (section  5.3.2)  as  well  as  the  stability  of  the  upstream  boundary  of  the  bubble  which 
supplies  the  wall- normal  momentum  needed  to  push  up  incoming  vortical  structures  (section 
5.3.2)  are  discussed.  Figure  5.15(d)  shows  the  turbulence  intensities  after  the  reattachment. 
The  SBL  case  displays  a  higher  degree  of  isotropy,  and  also  has  more  exaggerated  near  wall 
peaks  in  the  streamwise  and  spanwise  curves. 
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Figure  5.16  shows  the  turbulent  kinetic  energy  and  turbulent  kinetic  energy  production 
and  dissipation.  The  turbulent  kinetic  energy  is  a  compilation  of  the  rms  velocity  fluctu¬ 
ations  and  follows  the  similar  trend  of  the  turbulence  intensities  (figures  5.15(a)-(d)).  As 
the  detachment  point  is  approached  the  peaks  of  the  turbulent  kinetic  energy  and  turbulent 
kinetic  energy  production  and  dissipation  move  away  from  the  wall.  Both  the  production 
and  dissipation  decrease  in  magnitude,  while  the  turbulent  kinetic  energy  increases  (fig¬ 
ure  5.16(b)).  Inside  the  bubble,  the  turbulent  kinetic  energy  and  turbulent  kinetic  energy 
production  and  dissipation  become  smaller,  and  toward  the  reattachment  point,  the  vertical 
shifts  between  the  two  cases  become  negligible  (figure  5.16(c)). 

The  time-averaged  flow  field  is  similar  to  that  of  the  SBL  case.  There  exists  a  hump 
in  the  streamlines  throughout  the  boundary  condition  cycle,  which  remains  slightly  smaller 
than  the  separation  bubble  in  the  SBL  case  (figure  5.17).  The  contour  lines  shown  in  fig¬ 
ure  5.17  illustrate  the  maximum  and  minimum  heights  that  occur  during  the  boundary 
condition  cycle.  The  mean  streamwise  velocity  contours  above  u/U^  =  0.4  show  little  vari¬ 
ation  throughout  the  simulation.  Likewise  the  wall-normal  velocity  contours  in  this  upper 
region  maintain  a  symmetric  half-circle  shape  which  heaves  in  phase  with  the  syncopated 
Vtop{x,t).  Within  the  boundary  layer  the  mean  streamwise  velocity  contours  experience  a 
much  wider  variation  particularly  during  the  brief  Vtop,mm(x,t )  when  the  u/Uoo  =  0  col¬ 
lapses  onto  the  wall.  The  vertical  displacement  between  u/Uoo  =  0.0  and  0.4  maxima  and 
minima  is  significantly  different  (figure  5.18).  This  suggests  that  the  incipient  mechanism 
of  full  reattachment  is  a  thickening  of  the  shear-layer  which  blankets  the  separation  bubble. 
Streamwise  velocity  profiles  in  wall-coordinates  are  depicted  in  figure  5.19.  In  order  to  scale 
using  shear  stress  velocity,  uT,  during  separation,  Skote  &  Henningson  (2002)  ’s  definition  is 
employed: 


uT 


-V 


du 

dy 


y= o 


(5.5) 


Phases  5  and  9  are  chosen  because  they  represent  phases  in  which  the  maximum  and 
minimum  backflow  velocities  occur,  respectively.  Interestingly,  while  the  minimum  backflow 
velocity  occurs  at  the  phase  9  (approximately  —0.07  *  Ua c),  the  minimum  shear  stress 
occurs  during  phase  8  (figure  5.12(e)).  While  the  phase  in  which  flow  is  fully  attached  is 
easily  identifiably  using  =  1.0  and  Cf^max,  the  phase  in  which  the  separation  bubble 
is  fully  formed  is  much  harder  to  define.  In  fact,  the  maximum  minimum  skin- friction 
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coefficient,  minimum  backflow  velocity,  and  maximum  bubble  height  all  occur  in  consecutive 
phases,  in  the  order  listed. 

Vortex  identification  and  vorticity 

Vortex  identification  is  performed  using  the  A2  definition  of  Jeong  &  Hussain  (1995).  The 
large  scale  structures  present  are  consistently  elongated  in  the  streamwise  direction  (fig¬ 
ures  5.20(a)  and  (b)).  Their  distribution  before  and  after  the  separation  bubble,  or  hump, 
have  a  direct  influence  on  the  streaks  seen  in  wall  shear  stress  contours  (figures  5.20(c)  and 
(d)).  As  the  circulation  reestablishes  in  the  hump  following  following  full  reattachment, 
the  vortices  in  the  detachment  region  begin  to  bunch  up  like  a  rug.  This  stall  in  motion 
results  in  a  gap  in  structures  on  the  leeward  side  of  the  bubble  5.20(b).  Once  development 
of  the  bubble  is  complete,  the  vortices  advect  over  the  bubble  just  as  they  do  during  the 
SBL  case.  When  the  circulation  seizes  and  the  flow  reattaches  the  height  of  the  structures 
is  not  greatly  affected.  However,  downstream  of  the  hump  the  spanwise  spacing  of  large 
structures  increases,  and  subsequently,  Az,  the  average  spanwise  spacing  of  wall  shear  stress 
contour  streaks  (figure  5.20(c)).  Throughout  the  simulation  there  is  a  void  of  vortices  with 
the  bubble  and  hump,  despite  full  reattachment. 

The  orientation  and  inclination  of  the  vortices  indicates  that  the  magnitude  of  the 
streamwise  vorticity  dominates  the  vorticity  vector.  This  is  true  in  all  but  the  detachment 
region.  During  bubble  development  the  bunching  of  the  vortices  described  above  also  affects 
their  orientation,  and  thus  the  wall-normal  vorticity  is  actually  fed  by  the  structures  that 
were  oriented  streamwise  when  they  developed.  Instantaneous  contours  of  streamwise  vor¬ 
ticity  at  two  different  streamwise  stations  are  shown  in  figure  5.21.  Contrary  to  the  variation 
of  bubble  height  observed  based  on  mean  velocity  fields,  the  vorticity  field  in  the  leeward 
side  of  the  hump  actually  grows  in  height  upon  full  reattachment  (figures  5.21(a)  and  (b)). 
However,  downstream  of  the  bubble,  this  relationship  is  inverted  (figures  5.21(c)  and  (d)). 
This  is  best  explained  using  Taylor’s  Hypothesis  (Taylor,  1938).  The  heightening  of  the 
vorticity  distribution  seen  during  the  attached  phase(s)  advects  and  appears  downstream 
during  subsequent  phases.  This  also  indicates  that  the  region  of  decreased  momentum  which 
defines  the  hump,  begins  to  advect  during  full  attachment.  It  is  a  slow  process,  and  begins 
weakly,  so  this  motion  is  not  readily  apparent  from  mean  velocity  contours  or  large  scale 
structures. 

The  lifting  of  the  vortices  seen  in  figure  5.20  can  also  be  seen  in  sideviews  of  the  vorticity 
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distributions  (figure  5.22).  Na  &  Moin  describe  the  separation  bubble  in  the  SBL  case  as 
acting  as  a  streamlined  obstacle,  not  allowing  the  separation  zone  to  be  penetrated  by 
the  higher  momentum  flow.  In  figure  5.22(b)  one  can  see  that  the  seizure  of  circulation 
immediately  allows  the  hump  to  be  penetrated,  near  x/5*n  =  210.  Again,  only  weak  vortical 
structures  are  allowed  to  pass,  so  this  phenomenon  is  not  readily  apparent  when  viewing 
the  mean  contours  or  large  scale  structure  motion. 

5.4  Summary 

Direct  numerical  simulations  of  turbulent  boundary  layers  with  steady  and  unsteady  pres¬ 
sure  gradients  were  performed  to  improve  the  understanding  of  unsteady  separation  pro¬ 
cesses  of  turbulent  boundary  layers.  A  time  varying  blowing-suction  velocity  distribution 
along  the  upper  boundary  induces  an  unsteady  adverse  pressure  gradient  to  the  turbulent 
boundary  layer.  This  produces  a  separation  bubble  in  which  the  backflow  cycles  between 
full  separation  and  complete  seizure  of  circulation.  The  distinct  characteristics  of  unsteady 
separating  turbulent  boundary  layers  were  revealed  by  a  systematic  comparison  with  steady 
attached/separated  turbulent  boundary  layers. 

In  the  USBL  case,  the  process  of  creating  and  annihilating  circulation  is  observed  to 
be  hysteretic.  As  the  adverse  pressure  gradient  subsides  the  shear  layer  which  blankets 
the  separation  bubble  begins  to  thicken.  This  continues  until  the  flow  is  fully  attached. 
As  the  circulation  reappears  the  shear  layer  compresses  at  a  slower  rate.  As  a  result  the 
backflow  exists  in  a  flatter  region  and  relatively  high  negative  shear  stresses  are  observed. 
Throughout  the  simulation  the  displacement  thickness  in  the  region  of  the  bubble  remains 
large.  When  circulation  stops  the  bubble  does  not  appear  to  advect  downstream.  The 
hump  does  however  lose  some  of  its  strength,  and  allows  vortical  structures  to  penetrate  its 
upstream  face. 
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Figure  5.6.  RMS  velocity  fluctuations  in  the  APG  case  in  three  stream-wise  locations. 
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(a)  Velocity  profiles  upstream  of  the  separation 


(b)  velocity  profiles  downstream  of  the  reattachinent 


Figure  5.7.  Mean  streamwise  velocity  profiles  in  the  SBL  case  in  streamwise  locations  in  wall  units. 
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Figure  5.8.  RMS  velocity  fluctuations  normalized  with  the  freestream  velocity  in  the  SBL  case. 
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Figure  5.9.  Contours  of  the  instantaneous  skin-friction  coefficient.  Negative  values  are  dashed. 


Figure  5.10.  Amplification  factor  used  to  convert  Vtop{x)  to  Vt op(x,t)-  x-axis  indicates  the  phase 
number. 
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(a)  7u  distributions  in  phases  2-5. 


(b)  distributions  in  phases  5-8. 

Figure  5.11.  Spanwise-averaged  fraction  of  time  (7^)  that  flow  moves  downstream. 
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(a)  Phase  4 


(b)  Phase  5 


(c)  Phase  6 


(d)  Phase  7 


(e)  Phase  8 


(f)  Phase  9 

Figure  5.12.  Streamwise  distributions  of  phase  averaged  skin- friction  and  pressure  coefficients. 
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(b)  SBL 

Figure  5.13.  Streamwise  distributions  of  total  averaged  skin-friction  and  pressure  coefficients. 
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Figure  5.14.  Streamwise  velocity  profiles  in  wall  units.  —  USBL; - SBL. 
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(a)  x/5*n  =  120  (b)  x/5*n  =  170 


(c)  x/8*n  =  270  (d)  x/5*n  =  320 


Figure  5.15.  Turbulence  intensity  profiles  normalized  with  U^.  —  USBL; - SBL. 
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(a)  x/5*n  =  120  (b)  x/6*n  =  170 


(c)  x/ <5*„  =  270  (d)  x/SU  =  320 


Figure  5.16.  Comparisons  of  the  turbulent  kinetic  energy  and  turbulent  kinetic  energy  production 

u 3 

and  dissipation  normalized  with  100^,  .  —  USBL; - SBL. 
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Figure  5.17.  Contours  of  u/Uoo  =  0.4; —  SBL;  —  Phase  1;  —  Phase  6. 


Figure  5.18.  Contours  of  u/Uoo  =  0.4  and  u/Uoo  =  0.0;  Phase  1;  —  Phase  6. 
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Figure  5.19.  Mean  streamwise  velocity  normalized  with  S*r 
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*0. 


(c)  Contours  of  rwau  during  full  reattachment 


(d)  Contours  of  rwau  during  bubble  development 


Figure  5.20.  Instantaneous  vortex  fields  defined  by  A2  =  —0.1,  and  corresponding  wall  shear  stress 
distributions. 


54 


CHAPTER  5.  DNS  OF  UNSTEADY  SEPARATION 


(a)  Phase  1 


(b)  Phase  5 


O') 

3  * 


40  20 


(c)  Phase  2 


z/§: 


(d)  Phase  5 


Figure  5.21.  Instantaneous  vorticity  contours  at  x/S*n  =  260  ((a)  and  (b)),  and  x/5*n  =  320  ((c) 
and  (d)). 
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Figure  5.22.  Instantaneous  vorticity  distributions. 
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A  fully  conservative  finite  volume  method 
for  incompressible  Navier-Stokes 
equations  on  locally  refined  nested 
Cartesian  grids 

A.l  Motivation  and  objectives 

A  signifiant  portion  of  research  in  the  field  of  computational  fluid  dynamics  has  been  focused 
on  improving  solution  accuracies  and  reducing  computational  costs.  Despite  the  massive 
increase  to  our  computational  capacity  over  the  last  decade,  the  scale  at  which  we  require 
highly  accurate  large-scale  simulations  has  grown  at  a  significantly  faster  rate.  Large- 
scale  simulations  today  employ  hundreds  of  billions  cells,  impose  significant  computational 
costs,  and  require  massive  amounts  of  memory,  oftentimes  exceeding  the  capability  of  many 
mainframe  computers. 

One  of  the  most  common  ways  to  lower  computational  costs  is  to  lower  the  total  num¬ 
ber  of  grid  points  being  employed  throughout  the  domain.  The  unfortunate  tradeoff  to  this 
approach  is  a  significant  loss  in  simulation  fidelity.  To  combat  this,  some  have  employed 
higher-order  numerical  methods  to  improve  solution  accuracy  while  still  lowering  the  number 
of  computational  cells  (Peng  et  ai,  2003;  Shiau  et  ai,  1999;  Sau  et  ai,  2004).  This  approach 
has  been  shown  to  be  effective  when  compared  to  low-order  schemes  primarily  due  to  the 
cost  savings  associated  with  fewer  computational  cells.  However,  higher-order  methods  are 
oftentimes  very  difficult  to  implement  on  the  complex  grids  required  for  practical  applica¬ 
tions.  Furthermore,  higher-order  methods  tend  to  have  wider  discretization  stencils  making 
them  particularly  ill-suited  for  simulations  of  physical  phenomena  with  sharp  local  gradi¬ 
ents.  In  such  cases,  the  higher-order  schemes  effectively  serve  to  smooth  out  the  solution, 
yielding  a  solution  that  no  longer  accurately  represents  the  intended  physical  phenomena. 

A  more  practical  alternative  to  higher-order  methods  involves  effective  mesh  utilization. 
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Localizing  cells  in  areas  of  interest,  namely  regions  where  high  wave  numbers  are  present, 
can  reduce  the  total  number  of  computational  cells  required.  This  approach  is  especially 
important  in  flow  configurations  with  walls  and  obstacles,  as  the  resolution  must  be  sig¬ 
nificantly  high  near  these  objects  to  fully  capture  the  boundary  layer  physics.  In  simple 
configurations,  stretched  grids  generally  allow  for  effective  cell  utilization.  However,  as  flow 
configurations  become  more  complicated,  stretched  grids  become  less  effective  and  more 
time  and  effort  is  required  to  generate  grids  for  specific  computational  configurations. 

A  particularly  attractive  local  refinement  methodology  is  the  nested  Cartesian  grid  ap¬ 
proach  (Berger  &  Oliger,  1984).  In  this  case,  local  refinement  is  achieved  by  embedding  fine 
Cartesian  grids  within  coarse  Cartesian  grids.  The  attractive  properties  of  such  a  scheme  are 
that  Cartesian  grids  are  simple  to  create  and  manipulate  while  allowing  for  simple,  efficient, 
structured  algorithms  to  be  implemented.  Significant  prior  work  as  been  done  to  develop 
nested  Cartesian  grid  methodologies  (Berger  &  Oliger,  1984;  Khokhlov,  1998;  Durbin  & 
Iaccarino,  2002;  Pember  et  al. ,  1995;  Gerritsen  Sz  Olsson,  1998).  Unfortunately,  one  draw¬ 
back  to  using  nested  Cartesian  grids  is  that  simulations  are  generally  limited  to  Cartesian 
geometries.  Significant  work  has  been  done  over  the  last  three  decades  to  create  numeri¬ 
cal  algorithms,  such  as  immersed  boundary  methods,  that  can  handle  complex  geometries 
while  utilizing  Cartesian  meshes  (Peskin,  1972;  Kim  et  al.,  2001;  Tseng  &:  Ferziger,  2003). 
These  algorithms  have  been  combined  with  nested  Cartesian  grids  to  allow  for  simulations 
of  complex  geometries  (Roma  et  al.,  1999;  Ye  et  al,  1999;  Peng  et  al,  2010;  McCorquodale 
et  al,  2001). 

Even  though  locally  refined  grids  offer  many  computational  benefits,  numerical  errors 
introduced  at  mesh  refinement  interfaces  are  still  a  significant  concern.  Over  the  last  two 
decades,  significant  resources  have  been  invested  into  developing  methods  to  improve  the 
solution  accuracy  at  refinement  boundaries  (Almgren  et  al,  1998;  Minion,  1996;  Martin  & 
Colella,  2000;  Martin  et  al.,  2008;  McCorquodale  et  al.,  2001).  Interpolation-based  methods 
have  become  the  most  popular  answers  to  this  problem  because  they  can  be  tuned  to  control 
the  order  of  accuracy  at  mesh  interfaces  and  they  have  have  the  added  benefit  of  allowing 
for  efficient,  fully-structured  algorithms  to  be  implemented  for  each  individual  nested  block. 
Despite  this,  one  of  the  major  drawbacks  regarding  interpolated  schemes  is  that  they  are 
numerically  non-conservative  with  regards  to  kinetic  energy. 

Non-interpolated  algorithms,  on  the  other  hand,  allow  for  energy-conservative  numerics 
to  be  implemented  at  mesh  interfaces.  This  is  done  by  ensuring  that  face-centered  variables 
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are  defined  as  the  simple  arithmetic  mean  between  two  neighboring  cells’  values  (Mahesh 
et  al .,  2004).  Unfortunately,  non- interpolated  algorithms  on  nested  Cartesian  grids  suffer 
from  inherent  mesh  skewness  that  can  adversely  affect  the  solution  accuracy.  Prior  work 
on  reducing  the  adverse  effects  of  mesh  skewness  has  focused  on  correcting  the  numerical 
schemes  themselves  by  introducing  some  sort  of  skewness  correction  factor  (Mahesh  et  al, 
2004;  You  et  al.,  2008b;  Zwart,  2000).  However,  this  correction  factor  often  employs  some 
kind  of  weighted  interpolation  that  can  ultimately  destroys  a  scheme’s  ability  to  conserve 
energy. 

This  paper  will  introduce  a  new  hanging  node  scheme  designed  to  be  fully  conservative 
and  improve  on  prior  methods.  The  focus  of  this  method  is  to  geometrically  reduce  the 
inherent  mesh  skewness  found  in  nested  Cartesian  grids  by  virtually  slanting  cell  faces  at 
refinement  interfaces  to  align  cell  faces  with  the  normal  vector  connecting  neighboring  cell- 
centers.  This  implementation  will  allow  for  numerics  that  conserve  mass,  momentum,  and 
energy  while  retaining  the  benefits  of  an  overall  structured  implementation. 


A. 2  Computational  method 

A. 2.1  Fractional  step  method 

The  incompressible  momentum  and  continuity  equations  in  Cartesian  coordinates  are 

dui  duiUj  dp  d  dui 

dt  duj  dxi  U  dxj  dxj 

and 


(A.2) 


where  Ui,  p,  and  v  are  the  velocity,  pressure,  and  kinematic  viscosity.  For  incompressible 
flows,  the  density  is  considered  constant  and  is  absorbed  into  the  pressure  term.  The  prim¬ 
itive  variables,  velocity  and  pressure,  are  stored  in  a  collocated,  cell-centered  arrangement, 
with  an  independent  face  normal  velocity,  un  =  u  ■  n  stored  at  cell  faces. 

Time  integration  of  equations  (A.l)  and  (A.2)  is  performed  using  the  fractional-step 
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method  of  Kim  &  Moin  (1985).  The  semi-discretized  fraction-step  algorithm  is  as  follows: 


Ui  -  vJt 

dpn 

At 

dxj 

U*  —  Uj 

dpn 

=  A  t-£~, 

dxj 

d2pn+ 1 

1  u* 

dx2 

Ax  dxj  ’ 

<+1  -  < 


At 


duiUj 


dpn+1 

dxi 


+  za 


d2Ui 


dxY 


(A.3) 

(A.4) 

(A.5) 

(A.6) 


where  Uj  and  u*  are  intermediate  cell-center  velocities  and  the  convective  and  diffusive 
fluxes  are  integrated  using  the  third-order  accurate  Runge-Kutta  and  second-order  accurate 
Crank-Nicolson  schemes,  respectively.  All  spatial  discretizations  are  performed  using  non- 
weighted,  second-order  accurate  central  difference  scheme. 

Integrating  equation  (A.3)  over  a  computational  cell  volume,  AV ,  and  applying  the 
divergence  theorem  yields 


Uj  ~  Uj 

At 


AV 


Y2pnAS  -  ^2  UijaceunAS  +  ^2 
cs  cs  cs 


A  Uj 
Axn 


AS, 


(A.7) 


where  CS,  AS,  and  Ujjace  denote  the  control  surface,  the  surface  area  of  a  single  cell  face, 
and  the  interpolated  velocity  at  the  face. 


A. 2. 2  Conservation  principles 

Conservation  of  mass,  momentum,  and  energy  is  of  significant  importance  in  ensuring  that 
numerical  computations  capture  the  physical  phenomena  more  accurately.  For  compress¬ 
ible  flows,  discrete  energy  conservation  requires  that  the  convective  and  pressure  terms  are 
expressible  in  divergence  form.  The  following  derivation  is  available  in  the  literature  (Ma- 
hesh  et  al.,  2004)  and  is  repeated  here  for  clarity  as  to  the  implementation  in  the  current 
methodology. 

In  a  non-staggered  formulation,  this  is  complicated  by  the  fact  that  the  face-normal 
velocities  stored  at  the  face-centers  are  treated  as  independent  variables  in  relation  to  the 
Cartesian  velocities  components  stored  at  the  cell  centers.  We  can  see  this  by  looking  at 
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the  convective  term  for  a  passive  scalar  equation, 


^  +  ^-=0. 
dt  dxifiui 


(A.8) 


Integrating  equation  A.8  over  a  cell  and  applying  the  divergence  theorem  yields 


Vr 


d<t>c 


dt 


+ 


J2  Ha 


lAt  =  0, 


facesofcv 


(A.9) 


where  (f)cv  is  the  value  of  the  scalar,  Vcv  is  the  cell  volume,  and  Af  and  vn  denote  the  face 
area  and  the  face-normal  velocity,  respectively.  In  this  case,  since  the  incompressibility  con¬ 
dition  requires  that  =  0,  0  is  conserved  regardless  of  how  efface  is  computed.  However, 
conservation  of  <f>  does  not  imply  conservation  of  <^>2.  In  face,  <f> 2  is  only  conserved  if  (j)face 
is  computed  as  the  simple  arithmetic  mean  of  the  adjacent  cell-center  values, 

Vcv  ^ V  +  (fiface  =  ^  (H  +  ^nbr)  ■  (A.  10) 

volumes  boundary  faces 


Multiplying  equation  (A.9)  and  substituting  equation  (A.  10)  yields  the  discrete  equation 
for  (j)2, 

Vcv4>cv^^- +  ^2  (Hv  +  (frneigh)  VnAf  =  0,  (A. 11) 

facesofcv 

which  can  be  rewritten  and  simplified  as 

ir^k  +  ^T  V.  VnA,+\  Y.  ‘Pcv'PneighVnAf  =  0.  (A.12) 

facesofcv  facesofcv 

' - - - ' 

=0 

Summation  over  all  cells  in  the  domain  yields 

7/2 

^  ^  ^  ^  ^  '  4>cv4>neighvn-Af  =  0,  (A.  13) 

volumes  volumes  facesofcv 

in  which  the  contribution  of  all  interior  faces  cancel  out  and  result  in 

Vcv  +  ^2  4>cv(t)neighvnAf  =  0.  (A- 14) 

volumes  boundary  f  aces 
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As  long  as  4>neigh  is  also  defined  according  to  equation  (A.  10),  this  implies  that  both  4>  and 
cj)2  are  only  influenced  by  the  boundary  conditions  and  are,  in  face,  discretely  conserved. 
In  this  work,  discrete  energy  conservation  is  achieved  by  defining  the  interpolation  for  any 
face  variable,  </>jace,  such  that  it  is  the  simple  arithmetic  mean  of  the  current  cell  and  its 
neighboring  cell,  as  expressed  in  equation  (A.  10). 

A. 2. 3  Virtual-slanting  nested  Cartesian  grid  scheme 

Local  grid  refinement  is  implemented  by  introducing  logically  nested  Cartesian  blocks.  The 
benefit  of  this  approach  is  that  it  allows  for  simple,  computationally  efficient,  structured 
discretization  schemes  to  be  applied  with  a  lower  memory  requirement  compared  to  un¬ 
structured  algorithms.  Furthermore,  by  localizing  grid  refinement,  computational  cells  can 
be  localized  in  areas  of  interest  and  necessity,  ultimately  lowering  the  number  of  computa¬ 
tional  nodes  necessary  .  However,  these  nested  Cartesian  grids  carry  with  them  an  inherent 
skewness  which  can  adversely  affect  the  numerical  accuracy  of  the  discretization  schemes 
that  are  employed.  The  following  sections  will  describe  a  method  by  which  cells  located 
at  refinement  interfaces  are  virtually  slanted  to  reduce  this  inherent  skewness  and  improve 
both  the  numerical  accuracy  and  the  conservation  principles  of  a  numerical  scheme.  Fig¬ 
ure  A.l  shows  three  types  of  grids:  a  regular  Cartesian  grid,  a  non-slanted  hanging  node 
grid,  and  a  slanted  hanging  node  grid.  These  grids  will  be  used  to  illustrate  the  differences 
between  numerical  discretizations  on  each  grid  and  how  a  simple  geometric  relationships 
can  be  used  to  transform  grid  (figure  A.  1(b)  into  figure  A.  1(c)). 

For  the  sake  of  simplicity,  only  two  dimensions  will  be  considered  and  the  grid  spacings 
in  the  x-  and  y-directions  will  be  assumed  to  be  the  same,  namely  Ax  =  Ay.  The  grid 
resolution  of  uniform  cells  will  be  referred  to  as  Ax,  while  the  grid  resolutions  of  fine  and 
coarse  cells  will  be  referred  to  as  Axf  and  Axc,  respectively.  Furthermore,  the  coarse  and 
fine  grid  resolutions  are  assumed  to  differ  only  by  a  factor  of  two,  such  that  Axf  =  ^Axc. 

Regular  Cartesian  grids 

On  two-dimensional  regular  Cartesian  grids,  as  shown  in  figure  A.  1(a),  the  discrete  volume 
of  any  given  cell  is  easily  defined  as 


A  VCart  =  Ax  Ay  =  Ax2. 


(A.15) 
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Figure  A.l.  Three  example  grids:  (a)  a  regular  Cartesian  grid,  (b)  a  non-slanted  nested  Cartesian 
grid,  and  (c)  a  virtually-slanted  nested  Cartesian  grid.  It  can  be  seen  that  virtually-slanting  mesh 
(b)  aligns  the  face- normal  velocities,  un,  with  the  vectors  connecting  the  coarse  and  fine  cell-centers, 
thus  reducing  the  inherent  skewness  in  the  hanging-node  mesh. 
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Furthermore,  the  surface  area  connecting  any  two  cells,  AS,  is  defined  simply  as 


A  Scart  =  Ay  =  Ax 


(A.16) 


and  the  distance  between  the  cell-centers  projected  along  the  face-normal  vector,  Axn,  is 
equal  to  the  grid  spacing, 


A  xn^cart  —  Ax. 


(A- 17) 


Non-slanted  nested  Cartesian  grids 

For  the  nested  Cartesian  grids,  we  will  focus  only  on  the  cells  adjacent  to  a  refinement 
interface,  as  all  other  cells  are  covered  under  the  regular  Cartesian  grid  discretizations. 
Furthermore,  the  non-slanted  methodology  described  here  does  not  employ  interpolation  to 
fill  the  hanging  node  values.  Instead,  cells  at  mesh  interfaces  are  evaluated  directly,  much 
like  they  would  be  in  unstructured  implementations.  This  means  that  one  extra  face  must 
be  evaluated  for  coarse  cells,  while  the  number  of  operations  for  all  other  cells. 

On  non-slanted  nested  Cartesian  grids,  figure  A.  1(b),  we  can  see  that  the  coarse  and 
fine  cell  volumes  are  easily  computed  as 


A  Vf,NS  =  AxfAyj  =  Ax  j2 


(A.18) 


and 


A  Pc, ns  —  AxcAyc  —  Axc  , 


(A.19) 


just  as  they  are  on  the  regular  Cartesian  grids.  However,  the  surface  area  connecting  any 
two  cells  across  a  refinement  interface  is  now  defined  by  the  fine  resolution, 


(A. 20) 


Furthermore,  since  the  cell-centers  are  now  farther  apart,  the  face-normal  distance  between 
the  cells  becomes 


AxntNs  =  -A  Xf  =  -Axc. 


(A.21) 
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We  can  easily  observe  that  this  length  is  not  computed  along  the  vector  connecting  the 
two  cell-centers,  introducing  a  skewness  in  the  mesh  that  can  adversely  affect  the  numerical 
accuracy  of  the  discretization  schemes. 

Virtually-slanted  nested  Cartesian  grids 

To  reduce  the  inherent  skewness  at  refinement  interfaces,  the  faces  can  be  virtually-slanted 
to  align  their  normal  vector  to  the  vector  connecting  the  fine  and  coarse  cell-centers,  as 
shown  in  figure  A.  1(c).  The  amount  a  surface  needs  to  be  slanted  is  simply  a  function  of 
the  cell’s  geometry.  In  the  case  where  Ax  =  Ay,  the  appropriate  virtual-slanting  requires 
that  node  D  be  shifted  \Axf  away  from  the  coarse  cell-center.  The  new  fine  and  coarse 
cell  volumes  become 

A Vf>SL  =  AxfAyf  -  ^AxfAyj  =  ^Ax/2  (A. 22) 

and 

7  7 

A  VC,SL  =  -AxcAyc  =  -Axc2,  (A.23) 

b  b 

respectively.  Furthermore,  both  the  surface  area  and  face-normal  distance  between  cells 
change  such  that 

A  SSl  =  J  (jA^y  +  (Ay/)2  =  -d-Ax/  =  ^j^Axc,  (A.24) 

and 

A  xn,SL  =  ^Ax/)2+Qa  yfy  =  -d-Ax/  =  ^Axc.  (A.25) 

Aside  from  modifying  the  geometry  of  faces  along  a  refinement  interface,  it  can  be  seen 
in  figure  A.  1(c)  that  the  virtually-slanting  scheme  also  modifies  the  surface  shared  by  cells 
A  and  C.  The  regular  Cartesian  surface  area  of  Axf  between  cells  A  and  C  now  becomes 

1  2 

A  Sac  =  Ax/  -  A  x  shift  =  Ax  -  -Ax  =  -Ax.  (A. 26) 

Again,  this  is  a  simple  change  in  the  geometric  coefficients  for  these  cells.  Coarse  cells  are 
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not  affected  by  this  as  the  shift  only  occurs  within  the  fine  cell  block. 

Comparing  equations  (A.  18)  -  (A. 24),  we  see  that  the  only  difference  between  the  non- 
slanted  and  virtually-slanted  scheme  is  a  simple  modification  of  the  geometric  coefficients. 
Furthermore,  the  virtual-slanting  scheme  does  not  introduce  any  additional  neighbor  con¬ 
nectivity.  Therefore,  the  same  matrix  solvers  used  for  the  non-slanted  approach  can  be  used 
for  the  slanted  scheme  without  increasing  the  computational  cost. 

A. 2. 4  Virtual  slanting  limitation 

The  amount  of  virtual  slanting  necessary  to  align  a  face  to  the  normal  vector  connecting  the 
cell-centers  varies  according  to  the  cell  geometry.  Cells  with  a  high  aspect  ratio  will  require 
more  slanting.  However,  if  the  aspect  ratio  is  too  high,  it  will  not  be  possible  to  slant  a 
face  appropriately  without  propagating  mesh  skewness  into  the  nested  Cartesian  blocks  or 
without  creating  additional  cell  connectivity.  Figure  A. 2  shows  a  case  in  which  a  face  is 
slanted  to  the  point  of  creating  triangular  fine  cells.  It  can  be  seen  that  slanting  the  faces 
any  further  will  affect  the  fine  cells  located  on  the  interior  of  a  nested  Cartesian  block. 

One  can  avoid  creating  such  complicated  geometries  by  limiting  the  extent  to  which 
faces  are  virtually  slanted.  Assuming  the  coarse  and  fine  resolutions  differ  by  a  factor  of 
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two,  this  geometric  limitation  is  defined  by 


Axj_ 

A  Xi 


<  Vs 


(A. 27) 


for  any  given  direction  i.  The  slant  angle  could  be  fixed  such  that  for  all  directions  i  where 
Axj/Axi  >  y/S,  the  face  may  only  slant  to  the  point  where  Axj/Axi  =  VS.  While  this 
avoids  the  creation  of  complicated  geometries,  the  slanted  face  would  no  longer  aligned  to 
the  normal  vector  connecting  the  coarse  and  fine  cell-centers  and  the  overall  effectiveness  of 
this  scheme  will  be  lessened.  However,  even  partial  virtual-slanting  should  serve  to  reduce 
the  skewness  inherent  in  nested  Cartesian  and  improve  the  mesh  quality. 


A. 2. 5  Extension  to  three-dimensional  configurations 

In  the  case  of  three  dimensions,  slanting  the  fine  cell  faces  such  that  they  are  aligned  to 
the  normal  vector  connecting  the  cell-centers  generates  new  cell  connectivity.  Figure  A. 3(a) 
shows  a  diagram  of  the  mesh  topology  resulting  from  three-dimensional  mesh  slanting.  The 
shaded  regions  in  figures  A. 3(a)  and  A. 3(c)  highlight  the  new  faces  generated  in  the  case 
of  three-dimensional  slanting.  This  additional  connectivity  increases  matrix  sparseness  and 
increasing  the  overall  computational  cost  of  the  simulation.  Moreover,  these  new  faces  would 
also  require  their  own  slanting  in  order  to  align  them  with  the  normal  vector  connecting 
diagonal  cell-centers,  increasing  the  implementation  complexity  even  further. 

In  order  to  avoid  these  complications,  an  approximated  slanting  is  introduced  in  the  case 
of  three  dimensions.  Rather  than  virtually  slanting  each  face  to  perfectly  align  it  to  the 
normal  vector  connecting  the  cell-centers,  only  the  node  at  which  the  four  fine  cells  meet 
is  shifted.  This  creates  two  sub-faces  on  each  fine  cell,  each  of  which  is  perfectly  aligned 
to  one  of  the  two  components  of  the  normal  vector  connecting  the  two  cell-centers.  The 
dashed  lines  in  figure  A. 3(b)  illustrate  where  the  two  sub-faces  meet  on  each  fine  cell.  It  is 
interesting  to  note  that  the  creation  of  these  sub-faces  does  not  add  to  the  computational 
cost  of  the  algorithm  as  each  sub  face  does  not  need  to  be  evaluated  separately  for  every 
time  step.  The  reason  for  this  is  because  each  discretized  evaluation  calls  only  for  the 
projection  of  the  face  area  onto  the  normal  direction  under  consideration.  Since  no  two 
sub-faces  overlap,  the  projected  face  area  for  each  face  is  simply  the  summation  of  the 
projected  areas  of  each  sub-face  onto  the  normal  direction  being  considered. 
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Figure  A. 3.  Isometric  views  of  the  fine  cells  at  a  refinement  interface  are  shown  above  for  the 
(a)  fully-slanted  and  (b)  approximated  slanting  configurations.  A  side  view  of  each  each  method 
is  also  shown  in  (c)  and  (d),  respectfully,  for  clarity.  The  alternating  shaded  regions  in  (a)  and 
(c)  are  used  to  show  the  new  faces  generated  between  edge  cells  when  the  fully-slanted  approach  is 
employed.  These  new  faces  expose  the  coarse  cell  to  diagonal  cells,  requiring  additional  connectivity 
information  and  increasing  the  computational  cost  of  the  simulation,  (c)  and  (d)  clearly  illustrate 
that  no  new  faces  are  introduced  in  the  approximated  slanting  approach. 
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A. 3  Verification  and  validation 


A. 3.1  Unsteady  convection-diffusion  of  a  Gaussian  pulse 


Unsteady  convection-diffusion  of  a  Gaussian  pulse  is  evaluated  to  verify  the  accuracy  of  the 
present  virtual-slanting  scheme.  The  computational  domain  is  defined  on  [0,  2]2  with  the 
analytical  solution  to  the  problem  being 


4>(x,y,t) 


1  ^  /  (x  —  ut  —  0.5)2  (y-vt-  0.5)2 

At  +  1  6XP  (  v{At  +  1)  v{A t  +  1) 


(A. 28) 


The  analytical  solution  is  used  to  both  initialize  the  field  at  t  =  0  and  apply  Dirichlet 
boundary  conditions  at  the  edges  of  the  domain.  Cell  Reynolds  numbers  of  2  and  200  are 
considered  by  imposing  a  fixed  viscosity  of  v  =  0.01  and  setting  convection  velocities  of 
u  =  v  =  0.8  and  u  =  v  =  80,  respectively. 

The  computational  mesh  is  composed  of  two  different  refinement  levels,  as  shown  in 
figure  A. 4(a).  The  left  half  of  the  mesh  is  the  coarse  grid  with  a  resolution  of  Axc,  while 
the  right  half  of  the  mesh  is  the  fine  grid  with  a  resolution  of  A Xf  =  ^Axc.  Only  the 
resolution  in  the  x-direction  will  be  discussed,  as  Ax  =  Ay  is  assumed  for  all  results  in  this 
paper. 

Figure  A. 4(b)  shows  that  the  Gaussian  pulse  passes  through  the  mesh  refinement  inter¬ 
face  without  numerical  artifacts.  Figure  A. 5  shows  L2-norm  errors  with  respect  to  the  exact 
solution  for  the  present  method  at  Re  =  2  and  Re  =  200  and  various  levels  of  grid  refine¬ 
ment.  Grid  convergence  tests  were  performed  with  a  fixed  time-step  size  of  At  =  2.5  x  10-4 
for  Re  =  2  and  At  =  2.5  x  10-5  for  Re  =  200.  The  present  method  is  found  to  be  globally 
second-order  accurate  in  both  cases. 


A. 3. 2  2D  Taylor-Green  vortex 

Two-dimensional  Taylor-Green  vortices  are  simulated  to  both  identify  the  present  method’s 
accuracy  when  solving  the  incompressible  Navier-Stokes  equations  and  verify  its  energy 
conserving  properties.  The  computations  are  performed  on  a  [—1,  l]2  computational  domain 
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(a) 


Figure  A. 4.  Diagrams  of  (a)  the  computational  mesh  layout  employing  two  levels  of  refinement 
and  (b)  the  contours  of  the  convected  Gaussian  pulse  from  left  to  right  (at  t  =  0  and  t  =  1.0, 
respectively)  for  Re  =  2.  The  thick  black  line  along  the  centerline  of  (b)  denotes  the  refinement 
interface. 
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(a) 
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(b) 

Figure  A. 5.  The  above  plot  shows  the  L2-norm  errors  as  a  function  of  grid  spacing  for  (a)  Re  =  2 
and  (b)  Re  =  200.  Both  plots  show  that  the  virtually- slanting  scheme  is  globally  second-order 
accurate. 
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with  the  following  analytical  solutions  for  velocity  and  pressure: 

u  =  —  cos(nx)  sin(7ry)e_27r2*//Re,  (A. 29) 

v  =  sin(7rx)  cos(7ry)e_27r2t/Re,  (A. 30) 

p  =  — ^[cos(27rx)  +  cos(27ry)]e_47r2^i?e.  (A. 31) 

The  field  is  initialized  using  the  analytical  solutions  at  t  =  0  with  periodic  boundary  con¬ 
ditions  imposed  in  all  directions. 

The  computation  mesh  is  composed  of  quadrants  of  alternating  levels  of  refinement,  as 
shown  in  figure  A. 6(a).  This  arrangement  ensures  that  each  quadrant  is  surrounded  by  a 
refinement  interface,  even  across  periodic  boundaries. 

Grid  convergence  tests  are  performed  at  Re  =  100,  where  the  Reynolds  number  is  defined 
as  Re  =  UmL j v  with  Um  being  the  initial  maximum  velocity  and  L  being  the  vortex  length. 
A  fixed  time  step  of  At  =  0.001  is  used  to  ensure  that  CFL  <  1  for  all  cases. 

Figure  A. 7  shows  L2-norm  errors  for  various  levels  of  grid  refinement.  Again,  the  present 
method  is  found  to  be  globally  second-order  accurate  for  both  velocity  and  pressure. 

Figure  A. 8  shows  the  time  history  of  kinetic  energy  for  the  non-slanted  and  virtually- 
slanted  schemes.  The  kinetic  energy  is  evaluated  to  serve  as  an  indication  of  a  method’s 
ability  to  conserve  energy.  It  can  easily  be  seen  that  the  virtually-slanted  scheme  shows 
improved  kinetic  energy  conservation  properties  compared  to  the  non-slanted  scheme. 

A. 3. 3  Flow  in  a  lid-driven  cavity 

Flow  in  a  lid-driven  cavity  is  used  to  evaluate  the  virtually-slanted  scheme’s  ability  to  predict 
steady-state  solutions.  The  computational  flow  configuration  is  shown  in  figure  A. 9(a). 
Coarse  grid  spacing  of  Axc  =  A yc  =  0.0125  is  used  in  the  center  of  the  domain  with  local 
refinement  embedded  near  the  walls  using  A Xf  =  A yj  =  ^A xc  =  6.25  x  10-3 

The  computations  were  carried  out  at  a  Reynolds  number  of  Re  =  1000,  where  the 
Reynolds  number  is  defined  as  Re  =  UudL/v.  where  Uud  is  the  speed  of  the  top  lid  and  L 
is  the  length  of  each  side  of  the  square  cavity.  Contours  for  the  steady-state  solution  for 
the  lid-driven  cavity  are  shown  in  figure  A. 9(b)  and  centerline  velocity  profiles  of  u  and  v 
are  shown  in  figure  A.  10.  The  centerline  velocity  profiles  show  good  agreement  with  the 
results  obtained  using  a  uniform  mesh  with  resolution  Ax  =  Axf  =  6.25  x  10~3.  These 
results  demonstrate  that  the  present  method  is  able  to  achieve  comparable  results  to  uniform 
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(a) 


Figure  A. 6.  Diagrams  of  (a)  the  computation  mesh  layout  used  to  simulate  two-dimensional 
Taylor-Green  vortices  and  (b)  the  pressure  contours  of  the  initial  field.  The  thick  black  likes  in  (b) 
illustrate  the  refinement  interfaces. 
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(a) 


h 


(b) 

Figure  A.  7.  The  above  plots  show  that  the  virtually-slanted  scheme  is  globally  second-order 
accurate  for  both  (a)  velocity  and  (b)  pressure  at  Re  =  100. 
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Figure  A. 8.  The  kinetic  energy  time- history  shows  that  the  virtually-slanting  scheme  conserves 
energy  better  than  the  standard  non-slanted  hanging  node  scheme. 
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(a) 


(b) 

Figure  A. 9.  The  above  plots  illustrate  the  (a)  computational  grid  employed  and  (b)  the  u- velocity 
contours  for  the  steady-state  solution  at  Re  =  1000. 
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meshes  while  utilizing  fewer  computational  cells. 


A. 3. 4  Flow  over  a  square  cylinder 

Flow  over  a  square-cylinder  is  used  to  validate  the  present  scheme’s  ability  to  predict  un¬ 
steady  flow  phenomena.  The  computational  domain  is  defined  on  [0,  51Z4] 2 ,  where  D  is 
the  square-cylinder  diameter,  and  is  shown  in  figure  A. 11.  The  cylinder  is  located  in  the 
center  of  the  computational  domain,  ensuring  that  there  is  adequate  distance  from  both 
the  inlet  and  the  far-held  boundaries.  A  Dirichlet  boundary  condition  of  u  =  U  and  v  =  0 
is  imposed  at  the  inlet,  where  U  is  the  free  stream  velocity.  Free-slip  boundary  conditions 
are  imposed  at  the  far-held  boundaries  and  a  Neumann  boundary  condition  is  set  at  the 
outhow  boundary. 

Five  distinct  levels  of  local  rehnement  were  implemented,  as  shown  in  hgure  A. 11.  Again, 
each  level  of  rehnement  is  a  factor  of  two  different  from  it’s  neighboring  level,  with  the 
smallest  grid  size  being  A xj  =  7.8125  x  1CF3  near  the  wall. 

The  computations  were  carried  out  at  Re  =  100  using  a  time-step  size  determined  by 
the  CFL  =  1  criterion.  Figure  A. 12  shows  contours  of  the  instantaneous  vorticity,  while 
hgure  A.  13  shows  the  variation  of  the  lift  force  as  a  function  of  time.  Table  A.l  shows  the 
resulting  mean  drag,  RMS  lift,  and  Strouhal  number  (St)  for  the  present  scheme,  as  well  as 
from  a  variety  of  literature  results.  The  present  method  shows  good  agreement  with  data 
from  the  literature  for  the  Strouhal  number  and  RMS  lift.  The  mean  drag  appears  to  be 
under  predicted,  but  further  studies  are  being  performed  to  determine  whether  this  is  an 
issue  of  grid  resolution  or  a  limitation  of  the  virtually-slanted  scheme. 


St 

Mean  Drag 

RMS  Lift 

Virtual-Slanting  Scheme 

0.144 

1.36 

0.16 

Chung  &  S.-H.Kang  (2000) 

0.141 

1.46 

0.20 

Franke  et  al.  (1990) 

0.154 

1.61 

0.27 

Sohankar  et  al.  (1997) 

0.146 

1.48 

0.16 

Robichaux  et  al.  (1999) 

0.154 

1.53 

— 

Okajima  (1965) 

0.14 

1.45 

— 

Sohankar  et  al.  (1997) 

0.143 

— 

— 

Table  A.l.  Computational  results  for  the  Strouhal  number  (St),  mean  drag,  and  RMS  lift  at 
Re  =  100  for  the  virtually-slanting  scheme,  as  well  as  several  results  from  the  literature. 
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(b) 

Figure  A.  10.  Good  agreement  can  be  seen  between  the  virtually-slanting  scheme  and  the  uniform 
mesh  results  for  the  centerline  velocities  of  both  (a)  u  and  (b)  v.  The  uniform  mesh  resolution  was 
set  according  to  the  finest  resolution  in  the  locally  refined  case,  Ax  =  Ax/  =  6.25  x  10~3. 
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Figure  A.  11.  The  overall  computational  configuration  is  shown  above  in  (a),  with  a  closer  look  at 
the  refinement  interfaces  around  the  square-cylinder  shown  in  (b).  Five  distinct  levels  of  refinement 
are  shown,  with  the  highest  resolution  of  A Xf  =  7.8125  x  10-3  located  at  the  cylinder  wall. 
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Figure  A.  12.  Contours  for  the  instantaneous  z-vorticity  are  shown  in  (a),  with  a  closer  look  at  the 
contours  around  the  square-cylinder  presented  in  (b) . 
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Time 

Figure  A.  13.  The  above  plot  shows  the  smooth  variation  in  the  coefficient  of  lift  (Cl)  as  a  function 
of  time. 
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A. 3. 5  3D  Taylor-Green  vortex 

The  virtually-slanting  scheme  was  further  tested  using  the  three-dimensional  Taylor-Green 
vortex  case  shown  by  Shu  et  al.  (2005).  This  computational  configuration  is  known  to 
undergo  a  kinetic  energy  cascade,  whereby  the  kinetic  energy  in  the  system  drops  very 
rapidly.  The  length  of  time  required  for  this  kinetic  energy  cascade  to  take  place  is  a  strong 
indicator  of  the  amount  of  numerical  and  viscous  dissipation  in  a  discretization  scheme, 
where  a  longer  time-delay  indicates  a  lower  level  of  dissipation. 

Much  like  its  two-dimensional  counterpart,  the  computational  mesh  for  the  three-dimensional 
Taylor-Green  vortex  case  is  composed  of  octants  of  alternating  levels  of  refinement  and  is 
shown  in  figure  A.  14(a).  Again,  this  arrangement  ensures  that  each  octant  is  surrounded 
by  a  refinement  interface,  even  across  periodic  boundaries. 

The  computational  domain  is  defined  on  [0, 27r]3  with  periodic  boundary  conditions 
imposed  in  all  directions.  The  flow  field  is  initialized  at  time  t  =  0  using 


u(x,y,z )  =  sin(fex)  cos(fey)  cos(fe^),  (A. 32) 

v(x,y,z )  =  —  cos(kx)  sin(ky)  cos(kz),  (A. 33) 

w(x,y,z)  =  0,  (A. 34) 

P(x,y,z)  =  ^-[(cos(2z)  +  2)(cos(2x)+cos(2y))-2],  (A. 35) 


where  the  wavenumber  k  =  2ir/\  =  1.  Even  though  the  initial  conditions  impose  a  two- 
dimensional  flow  field,  the  existence  of  a  pressure  gradient  in  the  ^-direction  (shown  in 
figure  A.  14(b))  forces  the  creation  of  three-dimensional  flow  for  time  t  >  0. 

The  energy  cascade  described  previously  is  purely  a  numerical  phenomena,  therefore  a 
uniform  mesh  case  is  being  used  as  a  standard  for  comparison.  Figure  A.  15  shows  the  the 
time-history  of  kinetic  energy  for  a  uniform  mesh  case,  the  virtually-slanted  scheme,  and  an 
interpolated  scheme.  The  virtually-slanted  scheme  shows  good  agreement  with  the  results 
from  the  uniform  mesh  case,  both  of  which  appear  to  initiate  a  kinetic  energy  cascade  at 
very  similar  points  in  time.  Conversely,  the  interpolated  scheme  begins  rapidly  losing  energy 
much  sooner  than  either  of  the  previous  cases,  indicating  that  it  suffers  more  strongly  from 
numerical  and  viscous  dissipation.  These  results  strongly  suggest  that  the  virtually-slanted 
scheme  is  able  to  be  applied  to  three-dimensional  flows,  conserve  energy  quite  well,  and 
produce  results  comparable  to  those  found  on  uniformly  refined  meshes. 
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Figure  A.  14.  The  above  diagrams  depict  (a)  the  computation  mesh  used  to  simulate  three- 
dimensional  Taylor-Green  vortices  and  (b)  contours  of  the  initial  pressure  held. 
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Time 

Figure  A.  15.  The  above  plot  shows  the  time-history  of  the  kinetic  energy  for  three-dimensional 

Taylor-Green  vortices.  The  virtually- slanted  ( - red  line)  scheme  shows  good  agreement  with 

the  uniform  results  ( —  black  line),  while  the  kinetic  energy  cascade  begins  much  sooner  for  the 
interpolated  scheme  (-  •  -  blue  line). 
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A. 4  Summary 

A  second-order-accurate  finite-volume  method  for  kinetic  energy  conserving  simulations  of 
incompressible  flows  on  hanging  node  meshes  has  been  developed.  The  order  of  accuracy 
was  verified  to  be  globally  second  order  on  both  the  unsteady  convection-diffusion  equations 
and  the  incompressible  Navier-Stokes  equations.  The  present  scheme  showed  good  predic¬ 
tion  of  both  steady  and  unsteady  flow  phenomena  using  lid-driven  cavity  and  flow  over  a 
square  cylinder,  respectively.  Improved  energy  conservation  was  demonstrated  using  decay¬ 
ing  Taylor-Green  vortices  in  both  two  and  three  dimensions.  The  virtually-slanted  method 
showed  good  agreement  with  uniform  mesh  results  in  the  case  of  a  lid-driven  cavity,  and 
good  agreement  with  literature  data  in  the  case  of  flow  over  a  square  cylinder. 
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